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A  Strategy  For  Time  Dependent  Quantum  Mechanical 
Calculations  Using  A  Gaussian  Wave  Packet  Representation 
the  Wave  Function 


Shin-Ichi  Sawada ,  Robert  Heather,  Bret  Jackson 
and  Horia  Metiu 


Department  of  Chemistry 
University  of  California 
Santa  Barbara,  California  93106 


ABSTRACT 


We  develop  methodology  for  performing  time  dependent 
quantum  mechanical  calculations  by  representing  the  wave  function 
as  a  sum  of  Gaussian  wave  packets  (GWP),  each  characterized  by  a 
set  of  parameters  such  as  width,  position,  momentum  and  phase. 

The  problem  of  computing  the  time  evolution  of  the  wave  function 
is  thus  reduced  to  that  of  finding  the  time  evolution  of  the 
parameters  in  the  Gaussians.  This  parameteT  motion  is  determined 
by  minimizing  the  error  made  by  replacing  the  exact  wave  function 
in  the  time  dependent  Schroe'dinger  equation  with  its  Gaussian 
representation  approxlmant.  This  leads  to  first  order 
differential  equations  for  the  time  dependence  of  the  parameters, 
and  those  describing  the  packet  position  and  the  momentum  of  each 
packet  have  some  resemblance  with  the  classical  equations  of 
motion.  The . paper  studies  numerically  the  strategy  needed  to 
achieve  the  best  GWP  representation  of  time  dependent  processes. 
The  issues  discussed  are:  the  representation  of  the  initial  wave 
function:  the  numerical  stability  and  the  solution  of  the 
differential  equations  giving  the  evolution  of  the  parameters; 
and  the  analysis  of  the  final  wave  function.  Extensive 
comparisons  are  made  with  an  approximate  method  which  assumes 
that  the  Gaussians  are  independent  and  their  width  is  smaller 
than  the  length  scale  over  which  the  potential  changes.  This 
approximation  greatly  simplifies  the  calculations  and  has  the 
advantage  of  a  greater  resemblance  to  classical  mechanics,  thus 
being  more  intuitive.  We  find  however  that  its  range  of 
applications  is  limited  to  problems  involving  localized  degrees 
of  freedom  that  participate  in  the  dynamic  process  for  a  very 
short  time.  Finally  we  give  particular  attention  to  the  notion 
that  the  GWP  representation  of  the  wave  function  reduces  the 
dynamics  of  one  quantum  degree  of  freedom  to  that  of  a  set  of 
pseudo-particles  (each  represented  by  one  packet)  moving 
according  to  a  "pseudo-classical"  (i.e.  classical  like)  mechanics 
whose  "phase  space"  is  described  by  a  position  and  momentum  as 
well  as  a  complex  phase  and  width. 
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I.  IKTKODUCTION 

In  a  series  of  papers  Heller  developed  a  scheme  for 

computing  and  interpreting  time  dependent  quantum  mechanical 
processes  by  representing  the  wave  function  as  a  superposition  of 
Gauss.ian  wave  packets.  Since  each  packet  is  chara-cter ized  by 
several  parameters,  (the  position  and  the  momentum  of  the 
packet's  center,  a  complex  width  and  a  complex  phase)  the 
calculation  of  the  time  evolution  cf  the  wave  function  is  reduced 
to  that  of  the  time  evolution  of  these  parameters. 

In  his  applied  work  Heller  used  a  version  of  his  method 
(which  we  call  here  the  simplest  Heller  method  (SUM)  which  is 
based  on  two  simplifying  assumptions.  (1)  The  first  assumes  that 
if  we  must  represent  the  wave  function  by  a  sum  of  Gaussians,  we 
can  propagate  each  Gaussian  independently.  This  means  that  the 
equations  of  motion  for  the  parameters  of  a  Guassian  G^  are  not 
allowed  to  depend  on  the  parameters  of  another  Gaussian  Gg .  We 
call  this  the  independent;  Gaussians  approximation  or  IGA.  (2)  It 
is  furtheb  assumed  that  throughout  the  (i.e.,  collision  or  photon 
absorption)  process  the  width  of  each  Gaussian  is  smaller  than 
the  length  over  which  the  potential  changes.  This  allows  the 
use,  at  each  time  step,  of  a  second  order  Taylor  expansion  of  the 
potential  around  the  instantaneous  center  of  the  Gaussian.  We 
call  this  the  locally  harmonic  approximation  (LHA). 

SHM  was  used  successfully  by  Heller  to  analyse  a  variety  of 

time  dependent  processes  such  as  atom-diatomic  collision^,  photo- 
7  9 

dissociation  ,  photoabsorption  ,  Raman  scattering  and  atom 
diffraction  by  surfaces.^®  The  method  provides  accurate  results 
as  well  as  a  novel  and  beautiful  interpretation  of  quantum 
dynamics  in  terras  of  a  classical  language.  A  common  feature  of 
these  applications  is  that  they  all  deal  with  the  short  time 
dynamics  of  localized  quantum  degrees  of  freedom:  in  a  way  their 
success  reflects  mostly  Heller's  skill  in  identifying  important 
problems  that  fit  the  SHM  validity  conditions,  rather  than  the 


representation  oeyond  the  domain  in 
ihdoning  IGA  and  LHA.  This  is  done  in 
•k  and  requires  mostly  a  revision  of  th'e 
are  implemented. 


nsider  one  degree  < 
ssed  in  future  worl 
epresented  as  a  sut 


The  re.waining  sections  a.re  concerned  with  more  practical 
matters:  the  initial  choice  of  the  Gaussian  representation 

('Section  IV),  the  final  state  analysis  (Section  V)  and  the 
numerical  stability  of  the  M£M  equations. 

It  is  our  feeling  that  the  use  of  a  Gaussian  wave  packet 
representation  as  implemented  here,  is  likely  to  be  very  useful 
in  treating  quantitatively  problems  in  which  localized  quantum 
degrees  of  freedom  are  involved  in  dynamic  processes  of 
moderately  long  duration.  It  is  particularly  suited  to  problem? 
in  which  such  degrees  of  freedom  are  coupled  to  a  large  number  of 
classical  variables  whose  state  is  specified  only  statistically 
(e.g.  through  a  temperature)  since  the  coupling  of  classical  and 
quantum  degrees  of  freedom  presents,  in  this  framework,  no 
conceptual  difficulty. 
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II.  PROPAGATION 

II. 2  Ma-thenatical  Prel iraiHari es 

ri.2.A.  The  minimum  error  method  (MEM) 

We  are  concerned  with  equations  of  the  form 

iji  =  Aili  (II  .  1) 

where  i|i  is  an  unknown  vector  In  a  Hilbert  space  and  A  is  a  known 
time  independent  (this  restriction  is  not  necessajry)  operator. 

We  assume  that  we  know  a  physically  motivated  way  of  writing  ij)(t) 
in  the  form 

ij)(x:  t)  (t>(.x  :  (Xj^(t ) } )  (11.2) 

where  the  e.xplicit  functional  dependence  of  the  approximant  (j>  on 
the  parameters  ....  X^j  is  known  and  the  time  dependence  of  4) 

takes  place  exclusively  through  Xj^(t) . Xj^,(t).  Thus  we  can 

derive  the  time  evolution  of  4*  by  finding  the  parameter 
trajectories  X^(t)  which  satisfy 

3<t»/3t  =  X.  =  A<}>  .  ( II  .  3) 

oX  ^  1 

(Repeated  indices  are  summed  over).  Since  we  know  the  explicit 
dependence  of  <|) ,  A<{)  and  3<})/3X^  on  x  and  X,  we  can  use  (II. 3)  to 
develop  the  following  iteration  scheme.  We  assume  that  X^(t),  i 
=  1,  ....  N  are  known  and  use  (II. 3)  to  compute  X^,  i  =  l,  ....  N; 

then  we  can  determine  Xj,(t+T),  for  small  t,  from 

Xj^(t4-T)  =  X^(t)  +  X^T  O(T^)  .  (II. 4) 

and  repeat  the  procedure.  The  scheme  can  be  started  at  the 
initial  time  t=0  for  which  we  know  the  wave  function,  therefore 
the  values  of  Xj^(O),  i  =1,  ....  N. 
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The  equation  (11.3)  provides  us  with  an  infinite  num.ber  of 

equations  (one  for  each  value  of  x)  for  the  N  unknowns  Xj^(t).  To  . 

deal  with  this  situation  we  discretize  the  problem  by  using  Eq. 

(II. 3)  at  a  finite  number  of  points  x  ,  n=l,  ....  M ,  M  ^  N. 

n 

Thus,  we  have 

{3<}>(Xjj;{X})/3X^)X^  =  A4>(x^:{X})  (II. 5) 

If  we  denote  *  -3(|i ( x^ ) / 3X^  and  »  A(^(x^)  the  matrix  equation 

CX  =  B  is  a  set  of  M  linear  equations  with  N  unknown  and  M  ^  N. 
Such  equations  appear  in  the  "calculus  of  observations"^^ 
whenever  the  number  of  data  points  taken  by  overly  Industrious 
experimentalists  exceeds  by  far  the  number  of  unknowns  to  be 
determined.  A  customary,  but  not  unique,  way  to  get  the  "best" 
solution  is  to  minimize  the  quantity 

n 

with  respect  to  the  unknowns  X^.  The  weight  is  included  to 
allow  us  to  de-emphasize  the  role  played  by  the  less  reliable 
points  n,  or  to  enhance  the  influence  of  the  reliable  ones.  The 
extremum  (hopefully  a  minimum)  conditions  36/3X^=0  (for 
simplicity  we  assume  real  parameters)  lead  to  a  N  x  N  equation 

j  n 

+  + 
where  C  is  the  adjoint  matrix  of  C  and  C.  =  C  ..  This  equation 

in  ni  ^ 

has  a  solution  if  the  rank  of  the  matrix  '^n^ln*'n’j 

Since  D  is  the  Gram  matrix  of  C  the  rank  of  D  is  equal  to  the 
rank  of  C.  Thus,  the  parameters  Xjj^  can  be  determined  if  the  M  x 
N  matrix  =  3<{>  ( x^^  (X) ) /SX^,  has  one  N  x  N  minor  whose 

determinant  is  non-zero.  By  taking  the  continuous  limit  (the  x 
axis  is  divided  in  M  segments  of  equal  length  Ax^,  x^  is  taken  in 


the  middle  of  the  corresponding  segment  and  M  ->  e»)  and  using  W 


X  Ax  we  can  rewrite  (11.6)  as 
^n  n  ' 


/dxxtx)  (|^  X.  -  A^)*  (|^  X^  -  A<}>.) 


( II . 6b) 


and  (11.7)  as 


J‘dxx(x)  (|^)*  “  J‘dxx(x)  (f^)*  A()) 


{II. 7b) 


This  minimum  error  method  (MEM)  with  the  particular 
implementation  given  above  reduces  (when  we  take  x(x)  =  1)  to  the 


time  dependent  variational  principle  previously  used  in  quantum 
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mechanics.*"  The  chajige  of  the  point  of  view  introduced  by  the 
above  presentation  has  a  "liberating"  effect  since:  (a)  it  shows 
the  tremendous  richness  and  flexibility  resulting  from  the 
existence  of  a  large  number  of  legitimate  and  resonable 
definitions  for  the  "error"  S,  each  leading  to  different 
equations  for  the  propagation  of  Xj(t);  (b)  it  indicates  that 
this  is  a  mathematical  procedure  that  can  b.e  applied  to  the 
propagation  of  any  observable,  not  a  physical  principle  tied  to 
the  wave  function  and  the  time  dependent  Schrodinger  equation. 


Its  main  function  is  to  reduce  the  propagation  of  ij)(x,t)  in  the 

N 

Hilbert  space  to  the  computation  of  N  trajectories  ^j('t)  in  R  . 


II. 2. B.  A  Perturbation  Theory  Approach 


The  propagation  scheme  presented  above  seems  to  be  of  first 
order  in-the  time  step  t,  since  it  solves  Eq.  (II. 7b)  for  X  and 
then  uses  (II. 4)  to  f ind.  X^ ( t+T ) .  We  can  attempt  to  use  large 
time  steps  by  considering  that  (II. 7b)  is  a  first  order 
differential  equation  and  by  applying  the  Runge-Kutta  (RK) 
method.  However  if  Eq.  (II. 7b)  is  a  first  order  expression  in  T 
the  use  of  a  high  order  RK  procedure  would  be  incorrect,  for  the 
reasons  explained  below.  Let  us  consider  the  equation  X  =  f(t) 
O(t^)  and  compare  it  to  X  =  f(t).  The  RK  method  applied  to  X  = 
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f(t)  uses  the  expansion 

X(t+T)  =  X{t)  +  X(t)T  +  X(t)  t^/2  +  X’(t)  t®/3!  +  0(T^  ). 

( II .8a) 

«=  X{t)  +  f(t)T  +  f(t)T^/2  +  f  (t)T^/3!  +  O(t^)  {11.8b) 

However  If  the  equation  Is  X  =  f(t)  +  O(t^)  the  expansion  (II. 3b) 

2 

misses  the  third  order  term  tO(t  )  as  well  as  the  higher  order 

*  2  **  2 

terms  originating  from  0{t  )  and  0(t  ). 

To  check  whether  the  use  of  RK  method  to  solve  Eg.  (II. 7b) 
Is  legitimate  we  can  compute  (II. 8a)  by  perturbation  theory  and 
compare  it  to  the  equation  (II. 8b)  used  by  the  RK  method.  We 
find  that  the  two  procedures  coincide  only  when  a  certain 
definition  of  the  error  S  is  used. 


If  we  take  a  time  step  T,  causing  a  parameter  change 
6X(t)equal  to 

5X(t)  =  X(t)T  +  X  (t)T^/2  +  0(T®)  (II. 9). 

the  approximate  wave  function  changes  by 

64.(t)  =  <i>(t+T)  -  <j.(r)  =  (XT  ^  X  1^)  *  (Xt^^)2-0(t®) 

3X 

•  ■  (II. 10) 

The  same  change  can  be  written  as 

6(i.(t)  =  A<*j(t)T  +  I  A^<j)(t)T^  r  O(T^)  .  (II. 11) 

by  expanding  formally  (})(t  +  T)  =  exp{AT)  (}>(t).  Since  the  two 
expressions  must  coincide  order  by  order  we  have 


^1  = 


( II . 12a) 


and 
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8d>  •• 
3X  ^ 


A^<)> 


(II .12b) 


Th*5  corresponds  to  using  34)/3t  -  A4>  and  .3 ^(}>/3 =  A^4)  as  two 
independent  equations.  We  can  now  apply  the  minimum  error  method 
to  them. 


We  have  two  infinite  sets  of  equations  to  determine  2.\' 
unknowns  and  X^,.  Eq.  (II. 12a)  is  the  same  as  (II. 3)  (thus 
giving  the  false  impression  that  (II. 3)  is  valid  in  first  order 
only),  but  Eq.  (II.12b^  has  not  yet  been  used.  To  apply  MEM  to 
these  two  equations  we  define 

.  5  Xdx5<(x)[|^  X  -  A<|)]*[||  X  -  A<}>]  (11.13) 

and  use  3&^/3X  =  0  to  compute  X^,.  This  leads  to  Eq.  (II. 7b). 
Then  we  define 

Sg  2  J‘dxx(x)(||  X  -  <^(t)]*[||  X  -  <f(t)]  (11.14) 

with 

•Kt)  =  A^<|)  -  X^  (11.16) 

3X‘^ 

and  use  Z&^/dX  =  0  to  determine  X  (X  is  already  known  by  solving 
(II. 7b)).  This  leads  to 

«  * 

Xx(x)  1^  dx  X  =  Xdxx(x)  |~  $(x;t)  .  (11.17) 

Note  that  we  could  have  legitimately  defined  the  error  as  d  = 
and  used  3S/3X  =  0  and  3S/3X  =  0  to  generate  equations 
for  X  and  X.  The  equations  obtained  in  this  way  are  different 
from  (II. 7b)  and  (11.17),  In  particular,  they  do  not  give  for  X 
the  same  value  as  the  time  derivative  of  Eq.  (II. 7b). 
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'We  can  now  turn  to  our  original  question:  Is  X  computed  by 
taking  derivatives  of  Eq.  {II.7fa)  Identical  to  the  value  of  X 
given  by  Eq.  (11.17)?  A  straightforward  calculation  shows  that 
this  is  the  case.  Therefore  a  n-th  .order  RK  procedure  applied  to 
the  first  order  "variational"  equation  (II. 3)  i.s  equivalen^t  to 
the  use  of  a  n-th  order  perturbation  theory  within  MEM  and  is 
thus  wholly  justified.  This  is  a  pleasant  result  since  the  use 
of  the  existing  RK  programs,  which  compute  the  needed  derivatives 
internally,  can  save  a  large  amount  of  labor.  Note,  however, 
that  the  use  of  the  error  &  =  with  the  equations  3&/3X  =  0 

and  3a/3X  =  0  leads  to  equations  for  X  and  X  which  are  not 
equivalent  to  the  RK  expansion  of  Eq .  (II. 7b).  This  is  true  for 
other  error  definitions  that  we  have  considered. 

II. 2. C  A  Global  Hinlatua  Error  Method 


The  applications  of  MEM  discussed  so  far  were  all  made  by 
using  errors  defined  locally  in  time.  Below  we  discuss  an 
extension  of  the  method  which  has  a  truly  variational  character 
since  it  determines  the  trajectories  X^(t)  which  minimizes  the 
error  fufictional 


T 

a  =■  X  dtn(t) 
0 


A<j>]  ^1-  “  A,(j>] 

3Xj  J 


(11.18) 


The  local  method  var-ies  the  numerical  value  of  X^(t)  so  that  the 
error  a(t)  made  at  time  t ,  is  minimized.  In  (11.18)  the  whole 
curve  X^ ( t )  is  adjusted  to  give  a  minimum  value  to  a.  The  weight 
n(t)  has  been  incorporated  to  permit  us  to  emphasize  or 
deemphasize,  as  desired,  the  importance  of  some  of  the  points  on 
the  trajectory.  Taking  the  functional  derivative  3a/3Xj^(t')  and 
equating  it  to  zero  leads  to 


n(  t)Re<3<}>/3X^l  (  (34)/3Xj  )Xj-A<|) )  >  [6  { t-T) -6  (  t )  ] 


-2(dn/dt)Re<(3(})/3X  )  1  (  (34)/3,X.)X  -A(}>)> 

J  J 

-2n{t)Re<(34>/3X^)l[(3^4>/3Xi3>:j)X^Xj  +  X^ 

-2n(t)Re<A3(})/3X|^|  [  {34»/3Xj  )Xj-A4.]>  =  0 


(11.19) 


Tf\e  first  term  appears  because  we  have  not  specified  any 
constraints  for  the  variations  6Xj^(T)  and  6X^(0)  at  the  ends  of 
the  trajectory.  We  can  eliminate  it  by  taking  n{T)  =  n(0)  =  0. 
The  second  term  is  zero  if  n(t)  is  a  constant;  if  the  first  two 
terms  are  thus  eliminated  n(t)  disappears  from  the  equation. 


The  equation  obtained  above  is  rather  different  from  the 
preceeding  ones  and  there  are  no  theoretical  grounds  for 
rejecting  one  in  favor  of  the  other.  The  existence  of  so  many 
ways  of  generating  the  trajectories  in  the  parameter  space 
originates  from  the  fact  the  the  "best"  solution  of  an  infinite 
set  of  equations  having  a  finite  number  of  unknowns  is  not 
uniquely  defined. 


II. 3  PHYSICAL  APPLICATIONS  0?  MEM 

II. 3. A  The  propagation  of  the  wave  function. 

In  the  application  presented  below  MEM  is  identical  to  a 
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known  "variational  principle".  ’  Its  application  to  the 
propagation  of  Gaussian  wave  functions  has  been  briefly  discussed 
by  Heller. Therefore  we  present  it  here  with  a  minimum  of 
details  which  are  Indispensable  for  understanding  what  follows. 

We  consider  a  representation  of  the  wave  function  of  the 

form 

i|i(x:t)  =  Z  G^(x:{A^^  (11.21) 

A=  1 

where  the  parameters  A.  (t),  a  =  1,2,  ...,  C,  are  complex 

A«  A 

functions  of  time  and  T.  (t),  a  =  l,  2 . R.  are  real  functions 

Aa  A 

of  time,  and  C.  and  R.  are  integers.  The  equation  (11.21) 

A  A 

represents  the  wave  function  i})(x:t)  as  a  sum  of  localized 
"fragments"  G^  whose  explicit  dependence  on  the  parameters  A  and 
r  is  known.  For  a  variety  of  reasons,  well  summarized  in  Heller's 
papers, the  use  of  complex  Gaussians  for  G^  is  particularly 
advantageous.  Other  functions  of  the  form 

(Z  Q„(t)(x-p„(t))")G^(x:{A^^}.(r^^)) 
n  =  0 

where  Qj^('t)  and  Pjj(t*)  are  functions  whose  time  dependence  is  to 

be  determined  and  G.  is  a  complex  Gaussian,  have  similar 

A 

advantages  and  greater  flexibility. 

To  calculate  the  parameter  trajectories  (i.e.  the  time 

dependence  of  A.  and  T.  )  we  use  MEM  with  the  definition  (11.6b) 

Aa 

for  the  error  S  and  the  operator  A  =  (ih)  H,  where  H  is  the 
Hamiltonian. 
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The  use  of  the  weighting  function  x(x)  require  some  comment 
since  previous  applications  of  the  "variational  principle"  use 
x(x)=l.  In  m-ost  dynamic  problems  we  are  not  interested  in 
computing  the  wave  function,  but  its  projection  on  some  given, 
time  independent  state  |t(0)>.  For  example  in  computing  the 
total  absorption  cross  section  for  an  electronic  excitation  of  a 
molecule  by  light  we  must  propagate  in  time  the  nuclear  wave 
function  t(x:0)  =  <x|t(0)>  of  the  electronic  ground  state  on  the 
final  electronic  energy  surface  ( Franck-Condon  approximation  is 
Implied),  to  obtain  the  wave  function  T(x:t)  =  <x|T(t)>:  the 
Fourier  transform  of  <T(t)|T(0)>  with  respect  to  t  gives  the 
total  absorption  cross  section.  One  can  show  in  fact  that  such 
quantities  are  generalizations  of  the  one  parti''le  Green's 
functions  used  in  many  body  theory  from  the  case  of  a  quasi- 
particle  excitation  to  that  of  a  transition  from  one  many-body 
state  to  another.  If  our  Intent  is  to  compute  such  overlaps  we 

might  as  well  weight  the  error  &  accordingly  by  taking  x(x)  =* 

* 

t(x:0)  t(x:0).  Thus  rte  determine  and  so  that  the  wave 

functj.on  i))(x;t),  given  by  Eq.  (11.21),  has  minimum  error  for 
those  values  of  x  where  x(x),  hence  t(x:0)  is  not  zero.  If 
t(x:0)  is  very  localized  this  procedure  should  allow  us  to  fit 
i}i(x:t)  with  fewer  fragments  G^(  x ;  { A)  ,  (T)  )  than  in  the  case  when 
we  try  to  fit  the  wave  function  in  the  whole  space. 

Applying-  ME.'i  to  the  approximant  defined  by  (11.21)  and  the 
error  defined  by  (II. 6b)  gives 


with 


^Bp:A^  ®Bp;Aa^Aa  °Bp  "  ° 


{ReC„,  ^  +  Re(B^  „.A,  )  +  Re  E„.  =  0 

BbjAa  Aa  A|i|8b  Ap  Bb 


(II . 25a) 


( II . 25b) 


Previous  work  treated  all  parameters  as  if  they  were  complex  thus 
generating  one  unneeded  equation  for  each  real  parameter.  In  all 
the  cases  that  we  are  aware  of  this  does  not  lead  to  errors  or 
serious  complications  since  the  superfluous  equations  can  be 
eliminated  by  inspection;  they  are  linear  combinations  of  the 
other  equations.  For  more  complicated  representations  of  ili(x:t) 
it  is  easier  to  use  the  procedure  described  here,  which  gives 
only  the  necessary  equations,  thus  avoiding  the  extra  work  needed 
to  carry  out  the  elimination  mentioned  above. 


II. 3. B  The  propagation  of  various  observables 

Since  the  wave  function  contains  all  the  information  we  can 

possibly  want  to  know. about  the  system,  it  contains  superfluous 

information  whenever  we  are  interested  in  .a  small  number  of 

observables.  Assuming  that  there  is  some  proportionality  between 

the  amount  of  information  wanted  and  the  effort  required  to  get 
1 3 

it  we  might  hope  to  save  labor  by  determining  the  parameter 
trajectories  that  give  the  best  fit  to  the  observables  of 
interest  only,  rather  than  by  fitting  the  whole  wave  function. 

In  the  case  already  mentioned,  when  we  want  the  overlap  of  <}’(x:t) 
with  a  localized  function  t(x:0).,  we  can  hope  to  need  fewer 
"pieces"  (e.g.  Gaussians)  if  we  determine  i|i  ( x ;  t )  only  in  the 
spatial  region  where  t(x:0)  is  large.  Similarly,  if  the 
variation  of  the  wave  function  4>(x:t)  with  x,  at  a  fixed  post- 
collision  value  of  t,  has  a  broad  hump  on  a  length  scale  L  with 
small  wiggles  on  the  scale  1  superimposed  on  it,  then  a  matrix 
element  with  a  planar  wave  function  of  momentum  of  order  h2Tr/l  is 
totally  indifferent  to  the  existence  of  the  hump:  it  is  however 
very  sensitive  to  the  details  of  the  wiggles.  Therefore  a 
calculation  that  gets  the  wiggles  right  and  misses  the  hump  is 
quite  satisfactory.  Again,  one  can  hope  that  such  diminished 
demands  on  the  quality  of  the  wave  function  requires  less  work 
(i.e.  fewer  Gaussians)  than  the  case  when  we  attempt  to  fit 
4)  (  X  :  t )  with  wiggles,  humps  and  whatever  else. 

Since  MEM  is  a  method  of  solving  differential  equations, 
rather  than  a  variational  principle  specifically  tied  to  ihSifi/St 
=  Hi}i ,  we  can  apply  it  to  generate  parameter  trajectories  that 
give  adequate  results  for  some  observable.  Several  examples, 
which  should  provide  ample  illustration  on  how  to  proceed  in 
general,  are  given  below.  For  simplicity  we  confine  ourselves  to 
the  case  of  one  approximant  (rather  than  a  sum,  as  in  Eq . 

(11.21))  and  several  real  parameters  The  generalization  to 
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the  case  (11.21)  is  straightforward. 

.  (o:)  The  use  of  the  transition  probability  to  determine  the 
parameter -  trajectories. 

Let  us  assume  that  at  pbst-collision  times  we  are 
interested  in  the  probability  that  the  systems  described  by  >}>'(t) 
is  in  a  continuum  state  j k> .  Taking  the  time  derivative  of  the 
probability 

'P(k:t)  s  l<k|«Kt)>|^  ,  (11.26) 

using  Eq .  (II. 2)  to  approximate  <(>(t),  and  the  Schrodinger 
equation  to  compute  ■  we  obtain 


P(k:t)  =  (3P/3X^)Xj,  =  (2/h)  Ira  <k|H4)><4>lk>  s  f(k:t) 

(11.27) 

Since  this  must  be  satisfied  for  all  values  of  k  (spanning  the 
continuum)  we  have  again  an  infinite  number  of  equations  and  N 
unknowns  X^^.  MEM  determines  the  unknowns  by  generating  an  N.xN 
equation  for  them.  This  is  obtained  by  minimizing 


S  =  /dk  x(k)  {(3P/3X^)X^  -  f(k:t)} 


(11.28) 


with  respect  to  X^.  The  result  is 


(J-dk  x(k)  1^  |^)^j  =  /dkx(k)(|^)  f(k:t)  . 


(11.29) 


Since  we  know  the  functional  dependence  of  <}>  on  ^|(t)  we 

can  compute  the  matrix  elements  appearing  in  (11.29)  whenever  we 

know  the  values  of  all  the  X.  at  t.  This  allows  us  to  determine 

1 

Xj^(t)  and  X  ( t  +  T  )i:X^  ( t ) -i-X^  ( t )  T  ,  providing  us  wit-h  ah  iteration 
scheme  to  get  X^(t)  at  all  subsequent  times.  If  <{)  is  a  Gaussian 
and  ik>  a  planar  wave  then  <(|)!k>  is  a  Gaussian  and  <kiH<j)> 
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contains  moments  of  a  Guassian  (from  the  kinetic  energy  operator) 
and  integrals  of  the  form  /dx  exp  ( -ikx )  ( x-x  ^ { x  )  (|>  (  x  )' .  The 
latter  can  be  performed  analytically  if  V(x)  is-fitted  by 
Gaussians,  exponentials,  polynomials  or  any  combination  of  them. 

The  same  procedure  can  be  applied  when  we  are  interested  in 
2 

P^{t)  =  |<n|i}j(t)>l  ,  where  n  is  a  discrete  state.  The  equation 

of  motion  is 

(3P^/3X^)X.  =  (2/h)  Im<n  |Hc|)><4>|n>  h  fj^(t)  (II. tO) 

If  the  number  of  wa-ve  functions  ln>  used  in  Equation  (11.30)  is 
larger  than  the  number  of  unknowns  (this  is  always  the  case  in 
a  Hilbert  space  of  infinite  dimension)  then  MEM  gives  ' 

3P„ 

tS  X„(3P„/3X^)(3P^/3X.)]Xj  =  2  ^  f„(t)  .  (11.31) 

where  Xjj  is  a  weighting  factor. 

Note  that  the  projection  on  discrete  basis  sets  to  give  the 

probability  P^  is  of  Interest  in  bound  state  dynamical  problems 

(e.g.  a  semi-classical  external  field  drives  the  system  into  a 

steady  state  which  is  a  superposition  of  the  eigenstates  of  the 

system).  In  collision  theory  we  need  probabilities  of  the  form 
2 

I  <k  I  <n  I  i}i  ( t )  >  I  =  P  ;,(t)  in  which  k  describes  the  relative 
translational  motion  and  n  the  internal  states  of  the  fragments. 
MEM  can  be  applied  to  such  situations  (to  compute  trajectories 
determined  for  the  best  fit  of  these  probabiliites )  with  no 
additional  conceptual  difficulty. 

(p)  The  use  of  expectation  values  ter -determine  the 
trajectories  . 

If  we  are  interested  in  knowing  the  expectation  value  of  an 
operator  0  at  time  t  we  can  use  it  to  generate  trajectory 


equations .  We  have 


dP 

|^<4,(t)|OH.(t)>.ZO  .  Z  f„(t)  (II. 32) 

n  U 

where  P  is  defined  by  (II. 2C)  and  f  by  (11.27)  (we  assume  here 
n  •  n 

that  the  discrete  basis  set  provided  by  the  eigenvectors  of  0  can 
describe  adequately  the  dynamics,  thus  only  0^^  =  <nlO|n> 
appears).  V?e  can  now  use  the  error 

&  =  Z  (11.33) 

n 

in  which  the  probability  equations  (11.30)  are  weighted  by  the 
matri.'c  elements  of  the  opera'.tor  0.  Thus  the  importance  of  the 
states  likely  to  contribute  more  to  the  mean  value  of  0  is 
emphasized  and  the  others  are  weighted  down  or  multiplied  by 
zero.  'By  equating  with  zero  derivatives  of  S  with  we  obtain 
(11.31)  with  the  Weight  -x  =  <n|0|n>. 
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HI.  Approximate  Propavation  Schemes 
III.l  Introductory  Remarks 

The  equations  (11.25)  can  in  principle  be  used  to  find  the 
time  evolution  of  the  liave  function  by  finding  the  parameter 
trajectories.  If  the  physics  of  the  problem  fbrces  us  to  use  too 
many  Gaussians  we  might  have  to  abandon  the  method  or  to  look  for 
some  simplified  propagation  schemes. 

To  see  how  rapidly  the  complexity  of  the  method  can 
escalate  let  us  consider  a  time  dependent  quantum  mechv.nical 
problem  involving  two  three-dimensional  variables  R  and  r.  We 
need  nine  complex  parameters  for  the  width  matrix  for  the 
variable  R  and  nine  for  r;  we  must  also  use  terms  of  the  form  (R- 
Rt)*C  •(r-r^)  to  permit  correlations  between  the  two  degrees  of 
freedom  and  this  requires  nine  complex  parameters.  Thus  the 
characterization  of  the  width  of  the  Gaussian  requires  27  complex 
parameters.  To  this  we  must  add  6  positions,  6  momenta  and  a 
complex  phase.  Thus  if  we  deal  with  six  correlated  degrees  of 
freedom  we  need  a  total  of  68  real  parameters  per  Gaussian.  For 
ten  Gaussians  we  must  solve  680  first  order  differential 
equations . 

Assuming  that  in  the  dumbest  possible  way  we  saturate  the 
space  with  Gaussians  and  are  willing  to  solve  700  (or  even  7000) 
equations,  the  method  could  still  be  used  since  all  the  labor 
required  to  carry  out  the  integration  to  obtain  the  parameter 
trajectories  is  thus  roughly  comparable  to  ‘that  needed  in 
molecular  dynamics;  seven  thousand  equations  corresponds  there  to 
2333  atoms,  which  is  fully  within  the  capability  of  present  day 
computers . 

The  disadvantage  of  such  a  brute  force  attack  is  the  loss 
of  the  simplicity  which  makes  the  Heller  method  so  appealing  in 
the  first  place.  It  is  not  surprising  therefore  that  most  of 
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Heller's  effort  was  directed  towards  simplfying  the  parameter 
equations  of  motion.  Such  simplifications  are  physically 
motivated  and  their  success  depends  on  the  problem  being 
addressed.  Nevertheless  some  of  their  features  are  such  that  can 
be  discussed  in  a  general  setting. 


All  numerical  calculations  carried  out  so  far  have  used  two 

approximations.  (1)  If  the  wave  function  was  constructed  as  a 

sura  of  Gaussians,  it  was  assumed  that  matrices  .  , 

BpiAp  BpiAa 

ReC_.  ,  and  D.  .  are  diagonal  in  the  indeces  B  and  A  which 
Bb  ;  Aa  A|i ;  Bb 

label  different  Gaussians.  This  approximation  decouples  the 
Gaussians  and  we  call  it  here  the  independent  Gaussian 
appifoximatlon  (IGA).  (2)  If  we  assume  that  each  independent 
Gaussian  is,  throughout  the  collision,  narrower  than  the  spatial 
range  over  which  the  potential  changes  appreciably,  we  can 
further  simplify  the  matrix  elements  since  thf'V  can  be  evaluated 
by  expanding  the  potential  in  power  series  around  the  center  of 
the  Gaussian  and  by  retaining  the  first  three  terms  of  the 
expansion.  That  is,  at  time  t  we  use 


V(r)  V('?^)M3V(‘r^)/3?^)  (r-r^)  +  (l/2)32v{r^)/3r,.^  (r-rj.)^ 

where  "r^  is  the  center  of  the  Gaussian.  In  what  follows  we  call 

this  the  local  harmonic  approxiniution  (LHA).  The  main  appeal  of 

this  approximation  is  that  the  mean  position  (i.e.  the  center  of 

the  paeJ^et)  and  the  mean  momentum  of  packet  move  according  tj 

1 4 

classical  mechanics.  As  shown  by  Ehrenfest  this  property  has 
nothing  to  do  with  the  use  of  Gaussians  for  G^('r,{X^}):  it  is 
valid  whenever  the  region  over  which  G^  is  non-zero  is  smaller 
than  the  spatial  range  over  which  the  potential  changes.  A 
further  appealing  feature  of  LHA  is  the  fact  that  the  phase  is 

essentially  the  classical  action  along  the  trajectory  followed  by 
the  center  of  the  packet,  which  is  in  agreement  with  the  eikonal 
approximation. 
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The  theory  resulting  after  making  these  two  approximations 

is  called  in  what  follows  the  simple  Heller  method  (SHM).  It  has 

been  successfully  applied  to  a  number  of  problems  chosen  so  that 

the  risk  of  SHM  break-down  wa.s  minimized.  A  beautiful  example  is 

the  absorption  coefficient  of  a  pho todissoc ia t ing  molecule.  The 

initial  state  is  bound  and  very  localized.  The  absorption  cross 

section  is  given  (essentially)  by  the  Fourier  transform  with 

respect  tO  time  of  the  overlap  between  the  Initial  wave  function 

and  the  time  dependent  wave  function  obtained  by  propagating  the 

initial  wave  function  on  the  upper  state.  If  the  fragments 

produced  by  photo-dissociation  separate  very  quickly  (i.e.  they 

are  on  a  strongly  repulsive  potential)  the  overlap  becomes  zero 

very  quickly.  Therefore  they  need  to  propagate  a  very  localized 

packet  for  a  very  short  time;  it  is  not  likely  that  it  will  have 

time  to  broaden  to  the  extent  that  will  cause  LHA  to  give 

substantial  errors.  By  using  SH.M  Heller  has  developed  a 

beautifully  clear  and  simple  picture  of  the  connection  between 

the  absorption  spectrum  and  the  classical  motion  on  the  upper 

state  on  which  the  dissociation  takes  place.  That  stimulated 

1 5 

equally  elegant  experiments.  Other  successful  calculations 
involved  harmonic  oscillators  for  which  -  as  we  show  later  -  SHM 
is  exact. 

1  6 

Recent  calculations  by  Skodje  and  Truhlar  and  by  Heather, 
1 7 

Jackson  and  Metiu  show,  however,  that  the  method  fails  to  give 
correct  values  for  the  time  evolution  of  the  states  of  the  Morse 
oscillator.  We  are  thus  led  to  examine  both  theoretically  and 
numerically  the  two  approximations  mentio.ned  above.  Our 
conclusion  is  that  they  are  justified  only  under  special 
c i rcums  tances . 


III. 2  The  Local  Harmonic  Approxination  (LHA) 

1 1 1 . 2 . A  The  definition  of  the  approklmatlon 

l¥  ri  -  1-  X  .  ^ 

We  consider  here  the  case  of  one  normalized  Gaussian,  to 
isolate  the  effects  of  LHA  from  those  of  the  neglect  of  the 
interaction  between  Gaussians.  Using  the  equations  (A. 3)  ~  (A. 7) 
of  Appendix  A  we  can  write  the  MEM  equation  for  the  case  of  the 
Gaussian  approxlmant 

iKxjt)  ::  G(x)  =  exp{  (i/h)  [a(x-R)^  +  P(x-R)  +  y]  }  (III. la) 
as 


M^(0£+2a^/m)  +  M^tr-PR-iWm  +  P^/2m]  +  7^  =  0 

M,(a+2a^/m)  +  [y-PR-lha/m  +  P^/2m3  +  V  »0  , 
6  0 

Refl  *  0  , 

Imn  .=«  0  , 

with 

n  =  M2[2a(P/m-R)  +  P]  +  =  0  . 

*=  <(x-R)"  G|G>  , 

and 

Vjj  =  <(x-R)"g|VG>  . 

The  equation  (Ill.le)  leads  to 
R  =  P/ra  , 

and  this  together  with  (III. Id)  gives 

•  ,  i_  <G|VIG>  <GI  (3V/3x)  |G> 

“  2  "  "  3R  <G!G>  "  "  <G|G> 


(III. lb) 

( III . IC) 

(III. Id) 
(Ill.le) 

(  III  .2) 
(III  .3) 

(III. 4) 

(111.5a) 

(111.5b) 


The  equations  (III. 5)  are  more  general  than  the  procedure 
used  here  for  their  derivation.  Since  R  and  P  are  the 


r: 


expectation  values  of  the  position  and  momentum  operators  for  a 

Guassian^state  the  equations  (III. 5)  also  follow  from  Ehrenfesf s 
theorem . 


Combining  (in.ib)  and  (III.lc)  we  can  write: 
[r-PR-iK(o:/m)  +  P^/2m]  =“  (M4VQ-M2V2)/ 

and 

a+2cx^/m  =  ( «2''o"''2  ^  ^  ^ ‘'‘4 • 


(III .5c) 


(Ill.Sd) 


The  LHA  assumes  that  at  any  time  t  we 
potential  by 

V{x)  ::  V{R)  +  (  3  V  (  R )  /  3  R  )  (  x-R  ( t )  ) 

+  (1/2) (3^V(R)/3r2) (x-R(t))^  . 


can  replace  the 


(III. 6) 


Using  this  expression  for  V(x)  in  the  equations  (III. 5)  leads  to 
the  LHA  equations: 


X 

1  V. 

R  =  P/m, 

•  ^ 

P  =  -3V(R)/3R 

■V 

,  - 

V 

r  -  PR  -  iha/m 

ct  +2o:^/m  =  -  ( 1/ 2  )  3  ( R  ) /3R  ^  , 


(111.7a) 
(III. 7b} 

(111.7c) 
( III .7d) 


In  what  follows  we  attempt  to  establish  the  limitations  of 
the  LHA  equations  (III. 7)  by  comparing  them  to  the  ME.M  equations 
(111.5)  from  a  physical  and  a  numerical  point  of  view. 


N 


-2.8  The  magnitude  of  the  .  rrnr  made  bv  T.HA. 

Clearly  the  expansion  (III. 6)  is  valid  only  if  V{x)  is 


^  ■ 
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practically  constant  as  x  varies  around  R  over  a  spatial  range 
equal  to  the  width  of  the  Gaussian.  A  more  precise  statement  can 
be  made  by  retaining  the  next  terms  in  the  potential  expansion 
and  requiring-  thgt  they  should  contribute  less  than  10"  to  the 
equation  of  motion.  Adding  a  third  order  term  to  Eq.  (III. 6)  and 
using  it  in  Eq.  (111.5b)  we  obtain 

^^MEM  "  ^LHA^^^LHA  = f 1 ( t ) ^ / 4 ] [ 3 ( R ( t ) ) / 3 R ( t ) ^ ] [ 3 V ( R ( t ) ) / 3 R ( t )  ^ 

+  0(1^)  ( II I . 8) 

Here  is  ^iven  by  Eq.  (III. 7b)  and  P,.nw  is  the  ME.M  value  of  P 

li  H  A  M  £  M 

when  the  third  order  term  is  included  in  the  potential  expansion. 

1/2 

The  length  l(t)  =•  [h/2Imo:(t)]  is  the  width  of  the  Gaussian. 

The  error  made  by  using  P,„a  instead  of  is  less  than  10^  if 

JjxIA  MbM 

(1^/4) (3^V/3R®)/(3V/3R)  i  0.1  (III. 9) 

We  have  found,  by  a  similar  analysis,  that  the  errors  in 
the  other  LHA  equations  are  smaller  than  10"  if  (III. 9)  is 
satisfied.  In  other  words,  the  LHA  equation  (111.7b)  is  the  one 
giving  the  largest' error . 

"Xx 

For  an  exponential  potential  V(x)  =  e  ‘  Eq.  (III. 9)  gives 

£0.4  ( 1 1 1  .  10) 

and  for  a  repulsive  Lennard- Jones  potentiai 

1^  13-14/R^  <  0.4  ,  (III. 11) 

where  R^  is  the  value  of  R  at  the  turning  point  (where  we  expect 

LHA  to  have  more  difficulties).  For  a  kinetic  energy  of  0.05  eV 

(thermal  for  He),  cr  =  4A,  e/k  -  54°K  and  V(x)  =  46[(cr/x)^^  - 
0 

(ct/x)  ]  we  find  (from  (III. 11))  that  LHA  is  satisfactory  (within 
10^)  if  1  £  0.08  A.  Roughly  the  same  result  is  obtained  from 
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(III. 10)  for  X  ^  cs  0 . 2A  (a  rapidly  varying  "hard  wall" 
potential).  Less  stringent  conditions  are  required  in  the 
smo'other  regions  of  the  potentials. 

In  our  calculations  of  scattering  of  He  from  solid 
1  8 

surfaces  we  find  that  1  exceeds  this  value  in  all  cases,  even 
though  we  have  varied  the  initial  width  (both  Re  oc  and  Im  a)  in 
an  attempt  to  obtain  narrow  packets  in  the  interaction  region. 

To  get  a  better  understanding  of  the  breakdown  of  LHA  we 
carried  out  several  calculations  in  which  the  Gaussian  wave 
packet  (GWP)  (III. la)  is  propagated  in  the  Morse  potential 

V(x)  =  D(l  +  exp[-2X(x-x^)]  -  2  exp [ -X ( x-x^ ) ] }  ,  (I  II. 12 a) 

In  Fig.  1  w.e  plot  l(t)  as  a  function  of  time,  for  a  normalized, 
initially  narrow  low  energy  wave  packet.  Since  the  Morse 
potential  is  the  sum  of  two  exponentials  (one  of  which  has  the 
length  scale  (2X)"^)  the  validity  condition  for  LHA  is  given  by 
Eq.  (III. 10)  (with  X  replaced  by  2X) .  This  leads  to  IX  =  0.31. 

We  expect  LHA  to  work  best  either  for  a  low  energy  GWP,  which 
samples  the  lower  part  of  the  potential  which  is  nearly  harmonic, 
or  for  packets  which  are  initially  sufficiently  narrow.  We  see 
that  more  than  half  of  the  time  XI ( t )  is  above  IX  =  0.31, 
indicating  that  the  conditions  for  the  validity  of  LHA  are  not 
fulfilled. 

It  is  important'  to  realize  that  in  order  to  be  a  useful 
approximation  LHA  must  be  uniformly  accurate;  that  is,  if  fg(t) 
and  f  (t)  are  the  exact  and  the  approximate  values  of  a  parameter 

2i 

we  must  have  / | f ^ ( t ) -f , ( t ) | dt <e  as  well  as  max 

Co  X 

\t  (t)-f  (t)|<s-  for  te[0,T].  Here  T  ,is  the  time  interval  over 
which  we  need  to  know  the  evolution  of  the  packet,  and  and 
are  small  numbers  specifying  the  error  we  are  willing  to 
tolerate.  The  reason  for  this  can  be  understood  by  considering 


the  trajectory  of  the  center  of  the  packet.  Let  us  assiime  that 
LHA  gives  us  the  incori-ect  force  only  for  te[t^,  t^+A],  This 
will  distort  the  trajectory  for  the  remainder  of  the  t'ime,  even 
though  the  force  is  computed  correctly  at  all  t  £  t^  +  A,  because 
the  values  of  R(t^+A)  and  P(t^+A)  are  erroneous,  and  therefore 
the  trajectory  will  stray  from  the  correct  path  at  t  >  t^  +  A. 

A  more  precise  test  of  LHA's  accuracy  is  made  in  Fig.  2 
where  we  plot  -  and  -  3  V  (  R  { t ) ) /3R  ( t )  ,  which  are  the  right 

hand  sides  of  the  MEM  and  LHA  equations  (111.5b)  and  (III. 7b), 
respectively,  giving  the  evolution  of  the  mean  momentum.  Thus  we 
are  comparing  the  expectation  value  of  the  operator  -  3V/3x  to 
t-he  classical  force;  if  LHA  works  these  two  quantities  must  be 
equal.  Again,  we  see  that  this  is  not  the  case. 

In  evaluating  the  LHA  accuracy  we  must  keep  in  mind  that 
the  trajectories  of  R  and  P  are  not  measurable  in  a  quantum 
experiment.  Normally  we  measure  the  projection  of  the  asymptotic 
wave  function  on  a  set  of  final  states.  It  is  conceivable  that 
such  projections  are  not  very  sensitive  to  errors  in  trajectories 
and  LHA  might  be  better  than  an  analysis  of  the  trajectory  might 
suggest.  On  the  other  hand  these  trajectories  are  used  to  give  a 
qualitative  description  of  the  dynamic  process  in  a  language  that 
is  close  to  classical  mechanics;  large  errors  in  the  trajectory 
would  lead  to  a  misleading  qualitative  representation  of 
dynamics . 


1 1 1 . 2 . C  A  comparison  between  the  pseudo-classical  mechanics 
generated  by  MEM  and  the  classical'mechanics  (given  by  LHA). 

The  MEM  Equations  (III.5a-b)  have  some  resemblance  to  the 
classical  equations  of  motion  for  the  coordinate  and  momentum; 
when  LHA  is  used  they  reduce  to  Hamilton’s  equations  with  the 
classical  potential  V(R(t)).  To  emphasize  both  the  fact  that  the 
MEM  equations  (III  5.a-b)  are  quantum  equations  for  the 
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expectation  values  of  the  position  and  momentum  operator, in  a 
Gaussian  state  and  the  fact  that  they  resemble  classical 
equations,  we  call  them  pseudo-classical  equations;  and  we  use 
the  term  pseudo-classical  mechanics  for  the  motion  of  R(t)  and 
P{t)  generated  by.  them.  For  a  single  Guassian  the  difference 
between  MEM  and  LHA  is  thus  equivalent  (as  far  as  R(t)  and  P(t) 
are  concerned)  to  the  difference  between  the  pseudo-classical  and 
the  classical  mechanics.  Since  these  trajectories  are  used  to 
interpret  quantum  dynamics  in  a  pictorial,  classical-like 
language,  it  is  instructive  to  examine  them  in  detail. 

III. 2. Cl  The  Potentials 


The  "pseudo-classical  potential"  v  =  <G 1 V ] G>/<G | G> 
appearing  in  (III. 5b)  can  be  written  as 

V  =  'S  dy  exp[-y^]  V[R ( t ) +1  ( t ) y ]  (111.13) 

-« 

1/2 

where  l(t)  =  (h/2Imo:(t)]  is  the  width  of  the  packet. 

Since  the  greatest  contribution  to  the  integral  comes  from 

the  values  of  y  between  zero  and  one  the  center  of  a  packet 

located  at  R(t)  is  acted  upon  by  the  values  of  the  classical 

potential  between  the  points  R(t)  and  R(t)  +  l(t),  "averaged" 

2 

with  the  Gaussian  distribution  exp(-y  ).  A  more  precise 
statement  can  be  made  for  the  exponential  potential  V(x)  =  exp[- 
Xx]  for  which 

V  =  exp{-X[R-Xl^/4] }  .  (III. 14) 

Thus,  for  this  particular  case,  the  pseudo-classical  potential 

acting  on  the  center  R(t)  of  the  packet  is  equal  to  the  classical 

2 

potential  at  the  point  R-Xl  /4.  The  packet  moves  as  if  it  is  a 


^8 
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"ball"  with  a  tima  dependent  "radius"  XI  /4;  its  center  interacts 
with  the  potential  before  it  reae.hes  the  interaction  region  of  the 
classical  potential;  and'it  turns  around  before  It  reaches  the 
classical  turning  point.  Note  that  the  "radius"  of  the  "ball" 
varies  in  time  and  depends  on  both  the  width  of  the  Gaussian  and 
the  rate  of  spatial  variation  of  the  potential  at  the  site  where 
the  packet  is  located. 

The  physical  origin  of  this  behavior  is  the  same  as  that  of 

the  Heisenberg  uncertainty  relation.  The  "radius"  of  the  packet 
2  2 

is  given  by  1  /2  =  <G|(x-R)  |G>/<G|G>  which  is  a  measure  of  the 

quantum  fluctuations  of  the  position  operator  in  the  Gaussian 
state.  To  bring  the  classical  and  the  pseudo-classical  poten¬ 
tials  into  agreement  we  must  have  l{t)  •>  0,  which  means  Imcr>«>.  In 

2 

this  case,  however,  the  .mean  kinetic  energy  <G  |  { -li  /2m) 

2  2 

3  /3x  )|6>/<G|G>  becomes  infinite  and  so  does  the  expectation 
value  of  the  energy  operator.  This  happens  because  the  length  1 
and  the  quantum  fluctuations  of  the  momentum  are  related  through 
the  Heisenberg  relation  (Ap*l  5  h/'J  2  with  a  minimum  uncertainty 
equality  when  Rea=0). 

The  pseudo  potential  v  =  <G | V | G>/ <G j G>  corresponding  to  the 
classical  Morse  potential  (III. 12a)  is 

V  =  D{ 1 +exp [ -2X( R-x  -Xl^/2)]  -2exp(-X(R-x  -Xl^/4)])  . 

0.0 

We  compare  v  and  V  in  Fig.  3  for  various  values  of  Ima  (i.e.  1 
(t))  sampled  from  values  that  occur  in  the  -MEM  calculations. 

Since  l(t)  varies  in  time  in  the  course  of  packet  propagat'ion,  v 
is  time  dependent.  The  drawings  in  Fig.  3  show  the  instantaneous 
values  of  the  pseudo-potential  for  various  values  of  1. 

The  values  of  the  potential  and  the  pseudo-potential 
energies  as  a  function  of  time  are  shown  in  Fig.  4.  ‘  These  are 
obtained  by  propagating  R  and  o:  according  to  the  MEM  and  LHA 


equations,  respectively.  Thus  v  depends  on  the  MEM  values 
R{t)  and  Intc:(t),  while  V  depends  on  the  LHA  (i.e.,  classical) 
values  of  R(t).  It  is  important  to  note  that  the  dependence  of 
the  pseudo-potential  .V  on  lmc:.(t)  has  no  classical  analog.  The 
appearance  of  Imc<(t)  reflects  the  quantum  fluctuations  of  the 
position  operator  which  makes  the  average  value  of  the  classical 
potential  different  from  the  classical  potential  at  the  average 
position  (i.e..  <GjV(.x)!G>  pi  V{<6j.x;G>)  where  .\  is  the  position 
operator).  This  reflects  Heisenberg's  uncertainty  principle.  If 
we  want  to  think  of  the  p’seudo-classicai  motion  in  classical 
terms  we  must  accept  the  fact  that  the  variables  P  and  R  are 
coupled  to  a  "classical  time  dependent  field  Ima(t)"  whose  time 
evolution  is  prescribed  by  the  ME.M  equations  (111.5). 

The  graphs  in  Fig.  4  show  that  the  time  evolution  of  the 
pseudo-classical  and  the  classical  potential  energies  v  and 
V(R(t)),  respectively,  is  rather  different.  A  detailed  analysis 
inclicates  that  they  differ  both  because ‘.ME.M  and  LHA  give 
different  results  for  the  time  evolution  of  R(t)  and  because 
l(t),  which  enters  in  v  but  not  in  V,  varies  gr.eatiy  in  time.  It 
is  interesting  to  note  several  of  the  effects  of  the  "field" 
Imo:(t)  on  v  which  make  its  time  evolution  very  different  from 
that  of  V:  v  does  not  vary  periodically  in  time;  the  point  where 
V  is  ma.ximum  is  not  the  turning  point  of  R(t):  the  point  where  v 
is  minimum  is  not  the  point  of  ma.ximum  kinetic  energy. 

III.2.C2  The  classical  and  the  oseudo-classical  energies. 

Further  difference  between  the  pseudo-classical  and  the 
classical  mechanics  can  be  seen  by  e.xamining  energy  conservation 
in  the  two  theories.  Because  the  pseudo-classical  equations 
resemble  the  classical  ones  we  can  apply  to  them  the  procedure 
used  in  classical  mechanics  to  derive  the  energy  conservation 
condition.  That  is  we  multiply  P  =  -3v/3R  with  P/m,  replace  P/m 
by  R  in  the  right  hand  side  and  rewrite  the  equation  in  terms  of 
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a  total  time  derivative.  The  result  is 


dt 


V  ( R  )  } 


h 


3  V 

3 1  ma 


d  I  met 
dt 


{III .16) 


2 

The  quantity  P  /2m  -r  v  =  a(t)  is  conserved  only  if  dlmct/dt  =  0 
(i.e.  frozen  Gaussians)  or  3v/3lmc:  =  0. 


In  general  the  pseudo-classical  energy  a(t)  is  not 
conserved.  There  are  several  useful  ways  of  stating  the  reason 
for  this.  Since  the. pseudo-classical  potential  v  depends  on  the 
"external  time  dependent  field"  Imo:(t),  the  systems  of  equation.s 
(III5a-b)  is  not  conservative;  the  "particle"  (i.e.  the 
trajectory)  exchanges  energy  with  the  field.  Another  insight  in 
the  behavior  of’S(t)  is  gained  by  examining  the  total  quantum 
energy  of  the  state  G: 

<G;HiG>/<G!G>  =  a(t)  -  h;alV2mIma  .  (III. 17) 

« 

Since  this  quantity  must  be  conserved,  the  ps eudo -c  1  as s  i  ca  I 

energy  S(t)  varies  in  time  to  compensate  for  the  time  evolution 
2 

of  h.a;  /2m  Imo:.  The  latter  quantity  is  equal  to  <G  (?- 
-  2 

<G|PlG>)  ;G>/2m<GiG>)  which  is  the  momentum  fluctuation  in  che 
state  ;G>  (P-  denotes  the  momentum  operator  and  <g;p;g>  =  ?(t)). 
The  presence  of  this  term  in  the  total  energy  is  a  purely  quantum 
effect  which  reflects  Heisenberg's  uncertainty  principle.  The 
term  is  very  large  when  the  packet  is  localized  in  the  coordinate 
representation  . 

The  time  evolution  of  <H>  =  <G .  H  |  G>/<G :  G>  for  .MEM  and  LHA 
is  plotted  in  Fig.  5.  together  with  S(t)  given  by  the  two 
theories.  'iv'e  see  that  <H>  is  conserved  in  MEM  but  not  in  LHA, 
which  conserves  the  classical  energy.  In  many  cases  it  is  useful 
to  monitor  <H>  in  LHA  calculations  since  its  change  with  time  is 
a  fairly  sensitive  indication  that  LHA  is  breaking  down. 


31 


I I I. 3  THE  INDEPENDENT  GAUS5IAN  APPROXIMATION  (IGA) 

I  1 1  •  3  .  A  The  D.escr  iution  of  the  Appi*oximatlon 

This  approxiraatioi.’  is  obtained  by  cancelling  in  the  equa¬ 
tions  (A3-A7),  (giving  the  time  evo.lution  of  the  parameters)  all 

the  integrals  of  the  form  /dx,(  x-R .  I*' (  x-Rq  )  "^  G.G.  =  M(An|Bra)  and 
^  A  o  A  i3 

J*dx(x-R.)*’  G.(x)  V(x)G_(x)  =  V(An|Bo)  in  which  A  differs  from  B. 

A  A  D 

In  the  compact  matrix  notation  of  the  equations  (A.9)-(A.12)  this 
amounts  to  retaining  only  the  diagonal  part  of  the  matrix  M  and 
eliminating  the  V(An|Bo)  terms  (Ap^B)  from  v.  This  decouples  the 
components  of  the  vector  X  and  leads,  after  a  little  algebra,  to 
the  equations  (III.5a-d)  for  each  GaJissian,  There  is  thus  no 
coupling  between  the  parameters  belonging  to  different  Gaussians. 

The  IGA  achieves  a  considerable  saving  of  both  programming 
labor  and  computer  tine.  Its  validity  is  however  suspect  on 
physical  grounds.  On  one  hand,  the  assumption  i)i(x;t)  = 

Z  G.(x:(X(t)})  requires  the  Gaussians  to  add  up  coherently  to  the 

A  ♦ 

correct  wave  function  at  all  times,  while  on  the  other  hand  IGA 
eliminates  all  the  matrix  elemehts  through  which  the  time 
dependent  Sch roe dinger  equation  forces  the  Gaussians  to 
influence  each  other.  Unless  very  special  circumstances  are  at 
work,  it  is  hard  to  believe  that  independent  Gaussians  can  act  in 
concert  to  construct  an  accurate  expression  for  ij). 

III.3.B  The  role  of  the  coupling  between  Gaussians.  in  the 
pseudo-classical  mechanics. 

To  understand  the  implications  and  the  consequences  of  IGA 
it  is  useful  to  examine  the  role  of  the  neglected  coupling  from 
the  point  of  view  of  the  "mechanics”  controlling  the  motion  of 
the  center  of  the  packets.  We  continue  to  call  this  a  pseudo- 
classical  mechanics  even  though  when  a  sum  of  Gaussians  is  used 
to  represent  iji  ( x ;  t )  the  resemblance  to  classical  mechanics  is 
diminished.  The  motion  of  one  real,  physical  particle  whose  wave 


32 


function  is  described  by  a  sun  of  N  Gau^sians  is  represented  in 
the  resulting  pseudo-classical  mechanics  by  the  trajectories  of  N 
"pseudo-particles",  tracing  the  motion  of  the  centers  of  the 
Gaussians.  In  the  MEM  equations  these  trajectories  are  coupled 
to  each  other  and  to  the  "external  time  dependent  fields"  cc^  and 
(i.e.  the  width  and  the  phase  of  each  Gaussian).  IGA 
eliminates  the  coupling  between  the  trajectories:  LHA  eliminates 
the  coupling  to  the  "external  fields". 

In  what  follows  we  examine  the  motion  of  the  coupled 

pseudo-particles  by  considering  two  Gaussians  only.  That  is,  we 

consider  the  wave  function  »Ji(x:t)  =  6.(x)  +  G-(x)  satisfying  the 

'  ”  °  2  2 

time  dependent  Schrodinger  equation  with  H  =  -(h  /2m)V  ■t-V(x),  and 
look  at  the  equations  for  ,  Rg  and  Pg .  The  latter  are 

given  by  the  third  and  the  fourth  rows  of  the  matrix  equation  ^  = 
(^M)  ^*7  (see  Eqs.  (A.10-A.13))  and  are 

^“a^^A^"  "  K  “  ^3  (III. 18) 

and 

2o:g(Pg/m  -  Rg)  +  Pg  =  (III. 19) 

The  right  hand  sides  of  these  equations  are  very  complicated 
complex  functions  (through  M(An[Bm)  and  V(AnlBO))  of  all  the 
parameters  of  the  two  Gaussians.  Taking  the  real  and  imaginary 
parts  of  Eq.  (III. 18)  we  can  solve  for  R,  and  P.  to  obtain 

Ra  =  P^/m  -  (2Ima^)'^  ImF^  (III. 20) 

and 

=  ReFg  -  (Reoc^/Ima^)  ImF^  (IIl’.21) 

These  equations  'take  simpler  forms  under  the  conditions 
discussed  below.  The  one-Gaussian  terms  in  ^  and  v  (i.e.  ,  th^e 


terras  of  the  form  XG^(x)*f(x)  G^(x)  dx,  with  f{x)  a  real  function 

such  as  ( x-R . ) ” ( x-R_ )  or  ^x-R,)"v(x),  are  real  and  the  two 
A  D  ^  A  ' 

Gaussian  terms,  i.e.,  JG  (x)  f(x)  G  (x)  dx  with  A  B ,  are 

A  D 

complex.  Therefore,  if  we  neglect  the  imaginary  parts  of  the  two 
Gaussian  terms  we  obtain  ImF2=0  which  leads  to  (from  Eq. 

(111.20))  classical  form.  Furthermore, 

from  (III. 21)  we  obtain  P  »  ReF.  where  F^  depends  on  the  para- 
meters  of  all  the  Gaussians.  This  is  very  different  from  the 

4 

classical  equation  of  motion  P^=-3V ( R^) /3R^  and  from  the  pseudo- 
classical  equation  for  one  Gaussian  <G  |  3V/3x  ]  G>/<G  |  G>  .  A 

further  simplification  can  be  obtained  by  setting  all  two 
Gaussian  integrals  equal  to  zero  (i.e.,  we  make  the  IGA).  In 
that  case  we  obtain  the  equations  (III.5a-d)  for  each  Gaussian. 

A  complete  decoupling  occurs  and  each  Gaussian  evolves 
independently  according  to  the  pseudo-classical  equations 
(III. 5).  Thus,  within  the  IGA  the  wave  function  can  be  described 
as  composed  of  the  coherent  sum  of  packets  whose  centers  are 
m'ovi.ng  independently  on  the  potential  surface  according  to  the 
pseudo-classical  equations  of  motion  under  the  influence  of  a 
force  given  by  the  time  dependent  mean  potentials  v^  s 
<G  .  1  V  I  G  .  >/<G  J  G  ,>  and  v_  h  <Go  |  V|  G.XGj,  |  Go>  •  The  widths  and  the 
phases  of  these  Gaussians  are  also  uncorrelated. 

As  we  have  already  mentioned,  the  pseudo-classical 

mechanics  generated  by  a  two  packet  wave  function  deals  with  two 

"pseudo-particles"  moving  on  twocoupled  trajectories.  Each 

pseudo-particle  has  its  own  potentia.1  or  Vg,  which  depends  on 

the  time  dependent  fields  Ima  and  Imo:  ;  besides,  the  pseudo- 

A  o 

particle  A  is  acted  upon  by  forces  neglected  by  IGA  which  depend 
on  Pg“P^.  non- 

classical  nature  of  such  forces  is  obvious. 

It  is  interesting  to  note  that  in  the  early  days  of  quantum 
mechanics  it  was  popular  to  represent  the  time  evolution  of  one 
particle  wave  functions  in  terms  of  the  flow  of  a  continuous 


distribution  of  classical  like  "particles"  endowed  with  well 
defined  trajectories  and  momenta  and  interacting  through  an 
effective  stress  tensor.^”®  The  multiple  trajectory  pseudto- 
classical  representation  proposed  here  is  in  many  ways  similar  to 
the  representations  proposed  in  these  early  works. 

The  use  of  multiple  Gaussian  wave  functions  can  be  easily 
justified  by  the  greater  flexibility  (i.e.  larger  number  of 
parameters  in  the  leaSt  square  fit)  of  the  basis  set,  which  gives 
hope  for  greater  accuracy.  There  are  however  many  important 
situations  when  the  use  of  multiple,  coupled  Gaussians  must  be 
used  even  if  a  crude  but  qualitatively  correct  description  of  the 
scattering  process  is  desired.  This  happens  in  multiple  channel 
problems  in  which  the  channels  are  not  overlapping  in  either  the 
coordinate  or  the  momentum  space.  One  example,  provided  by 
surface-atom  scattering,  is  the  case  when  one  channel  is  a 
particle  trapped  at  the  surface  and  the  other  is  a  particle  back- 
scattered  into  the  vacuum.  Another  exam-ple  is  provided  by  the 
curve  crossing  problems  in  which  an  atom  in  the  "ionic"  channel 
has  in  the  classical  limit  a  different  momentum  than  the  atom  in 
the  "neutral"  channel.  Such  events  cannot  be  described  -  even 
qualitatively  -  by  one  Gaussian  packet.  Therefore,  in  such 
situations  the  multi-pseudo-particle  description  of  the  dynamics 
is  the  only  reasonable  "classical  like"  picture  of  the  quantum 
process . 

IH.3.C  The  validity  conditions  for  IGA.  . 

Given  the  great  simplification  introduced  by  IGA  it  is 
important  to  have  a  clear  idea  under  what  circumstances  we  expect 
it  to  work.  We  discu's  ,  first  the  cas.e  when  IGA  is  used  together 
with  the  local  harmonic  approximation  (LHA)  (we  call  this  the' 
simple  Heller  method  (SHM))  and  show  that  if  LHA  is  made  then  the 
■wave  packets  become  decoupled  and  IGA  is  exact.  We  consider  this 
to  be  a  rather  striking  result  since  we  could  not  find 


intuitively  any  link  between  tbe  'two  approximations:  one  of  them 
has  to  do  with  the  relationship  between  the  width  of  each 
Individual  Gaussian  and  the  rate  of  change  of  the  potential  with 
x;  the  other  with  the  overlap  between  different  Gaussians.  Since 
we  believe  that  LHA  is  likely  to  fail  in  some  (or  many?) 
practical  cases  the  above  observation  is  not  of  much  pract-ica,! 
help.  It  does  how.ever  explain  why  Heller  was  so  successful  while 
using  the  IGA  method  in  problems  involving  harmonic  oscillators. 

If  LHA  is  not  made  one  can  show  that  the  Gaussians  might  become 
decoupled  when  the  packets  do  not  overlap,  or  when  they  have  very 
different  momenta,  or  when  their  phases  vary  extremely  rapidly  in 
time  . 

III. 3. Cl  The  simple  Heller  method  (SHM  =  LHA  +  IGA) 

There  is  numerical  evidence  that  harmonic  oscillators  have 
special  properties  with  respect  to  the  Gaussian  propagation 
method  discussed  here.  One  of  the  very  first  calculations 
carried  out  by  Heller  studied  the  excitation  of  a  harmonic 
oscillator  hit  by  an  at.im.  He  described  the  initial  oscillator 
wave  function  as  a  sum  of  Gaussians,  used  SHM  to  propagate  them 

and  obtained  satisfactory  results.  One  the  other  hand  both 
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Skodje  and  Truhlar  and  Heather,  Jackson  and  Metiu  have  shown 

that  SHM  or  IGA  gives  inaccurate  results  when  applied  to 

propagate  states  of  a  Morse  oscillator.  This  is  th’e  case  even 

for  low  energy  states  which  are  nearly  harmonic. 

In  order  to  understand  why  Heller’s  calculation  was  so 
successful  we  have  investigated  the  effect  of  LHA  on  the  coupling 
between  Gaussians.  We  have  found  that  LHA  decouples  the 
Gaussians  exactly.  As  a  corollary,  in  the  case  of  a  harmonic 
oscillator,  where  LHA  is  exact,  the  simple  Heller  model  (LHA  and 
IGA)  is  exact! 


To  show  how  this  is  proven  we  consider,  as  an  example,  the 


equation  (A. 5).  For  two  Gaussians  this  can  be  written  as: 

(P^/2ra  -  -  ih«^/m  +  y H(B1|A0)  +  (IM.22) 

[2o:^(P^/n-R^)  +  P^]  M(BllAl)  + 

[2ag(Pg/m  -  Rg)  +  P^^]  M(B2|B0)  + 

[ct^  +  2a^/m]  M(B1|A2)  + 

V{B1 I  AO  )  +  V(B1 1  BO)  =  0  . 

The  notations  M  and  V  have  been  specified  in  the  Appendix  A. 

The  discussion  proceeds  now  as  follows.  Let  us  assume  that 

the  two  Gaussians  move  according  to  the  simple  Heller  model 

(SHM):  that  is  the  quadratic  approximation  of  the  potential  is 

made  and  the  -Gaussians  are  assumed  to  be  independent.  This  means 

that  we  assume  the  SHM  equations  (Ill’.Sa)  and  (III.7b-d)  for  the 

parameters  R.,  P.,  7,  and  a,  and  R_,  P_,  v.  and  a„.  If  these 
AA  A  A  BB  B  B 

equations  are  inconsistent  with  (III. 22)  then  by  introducing  them 
in  (III. 22)  we  must  obtain  a  non-zero  result  whose  magnitude 
indicates  the  e.ttent  of  the  error  made  by  SHM.  Making  the 
substi  tuti.on  just  mentioned  le-ads  to 

Error  =  -  V(R^)  1  (  x-R^ )  G  ^G^dx  -  (  3  V/SR^)  J- (  X-Rg )  (  x-R^ )  G*  G^dx 

-(3V/3Rg)X(x-Rg)^G*Ggdx  -  ( 1 /2  )  (  3^  V/3R^)  f  (  X-Rg  )  (  x-R^ )  ^G  *G^  d:< 

^/(x-RglGgVG^dx  +  J'(x-Rg)GgVGgdx  =  0  .  (III. 23) 

Now  let  us  make  the  local  harmonic  approximation  to  evaluate  the 
integrals  present  in  the  error  expression.  If  we  expand  V(x)  in 
the  last  integral  in  powers  of  (x-Rg)  and  V(x)  in  the  integral 
before  the  last  in  powers  of  (x-R.),  we  find  that  the  error  is 


exactly  zero  (if  we  retain  only  the  quadratic  ten><  :  '! 
expansion)^.  Therefore  once  we  accept  LKA  the  i:u;ep..-..c;?-  . 
Gaussians  follows! 

III.3.C2  The  General  Case 


While  in  general  the  coupling  between  Gaussians  must  bu 
taken  into  account,  there  are  several  situations  in  which  it  can 
be  neglected,  even  if  LHA  is  not  made. 


(a)  The  most  obvious  one  is  when  the  Gaussians  do  not 
overlap.  This  can  happen  when  dealing  with  problems  in  which  the 
wave  function  tends  to  split  into  spatially  separated  pieces.  A 
trivial  example  is  the  low  energy  state  of  a  double  well 
potential.  A  more  i-nteresting  one  is  provided  by  atom  scattering 
from  a  moving  surface.  There  is  a  finite  probability  that  during 
the  collision  the  incident  particle  excites  phonons  and  is 
trapped  at  the  surface:  there  is  also  a  finite  probability  that 
the  particle  is  scattered  back  into  'the  vacuum.  Therefore  the 
atomic  wave  function  "splits"  into  a  component  bound  to  the 
surface  and  an  outgoing  free  particle  compone.nt.  If  th.®  wave 
function  is  approximated  by  two  Gaussians  they  will  best  mimic 
this  situation  if  one  of  them  is  trapped  at  t'ne  surface  and  t.he 
other  is  reflected.  E.xcept  for  the  early  times  during  the 
collision,  when  nothing  much  happens,  the  overlap  between  these 
Gaussians  should  be  fairly  poor  and  a  calculation  ignoring  th.e 
coupling  between  them  has  a  fair  chance  of  success. 

(p)  Another  interesting  situation  takes  place  when  the 

integrands  in  the  qua,ntities  .M(AniBm)  and  V{An'BO)  appearing  in 

the  equations  of  motion  (A. 9-12)  {or  the  Equation  (III. 22)  which 

is  one  particular  example)  oscillate  very  rapidly  around  zero. 

Since  all  such  terms  are  of  tiie  form  jdx  G^  Gg  f(x)  with  f(x)  a 

real  function  (of  the  form  (x-R.)”  (x-Rj™  or  (x-R)^V(x))  the 

Ad  a  ^ 

oscillatory  behavior  arises  from  the  phase  of  the  product  G.  G_ . 

A  D 

If  the  wave  length  of  this  oscillation  is  much  smaller  than  the 

* 

width  of  the  Gaussian  G.  6^  the  integral  is  practically  zero. 


One  term  in  the  phase  of  G,  G.  is  ( ( P_-P . ) /R] x ,  which  gives  the 

A  B  BA  • 

wave  vector  k  =  (P_  -  P.)/h.  Since  the  width  of  G.  G^  is  I  = 

221^2°^ 

1  g/ ( 1  1  g)  '  ,  the  integral  tends  to  zero  if  2Tr/k  <<  1.  This 

is  easy  to  understand  on  physical  grounds.  If  P_  and  P.  are  very 
»  BA 

different,  the  packets  G  and  G_  are  segments  of  planar  waves 

A  o 

having  very  different  wavelength.  As  is  well  known  such  waves 
are  poorly  coupled,  which  means  that  their  matri.x  elements  ,“S 

A 

f(x)  G„  dx  are  very  small, 
o 

(y)  Finally,  we  point  out  that  it  is  possible  that  two 

Gaussians  become  decoupled  if  their  two-Gaussian  integrals 

oscillate  rapidly  around  zero  with  time.  To  explain  this  we  can 

use  Eq.  (111.22)  as  an  e.x ample.  The  integrals  M{An  Bm)  and 

V(An|B0)  are  complex  and  therefore  have  the  form  a(t)e 

where  a(t)  and  4>(t)  are  real  functions  of  time.  The  structure  of 

these  integrals  is  such  that  they  will  have  the  same  phase  since 

that  is  determined  by  G.  G_  which  appears  in  all  integrands.  To 

A  o 

simplify  matters  consider  the  schematic  representation  of  Eq. 
(III. 22)  provided  by 

a(t)e"-'*’^^^X,  (t)  -  .M(t)X2(t)  =  V(t)  -  e"  ^ ^ ^  b  (  t )  (!n.24) 

Here  the  terms  with  the  phase  <|>(t)  aire  two-Gaussian  integrals, 

M(t)  and  V(t)  are  one-Gaussian  integrals  and  Xj(t)  and  .X2(t)  are 

combinations  such  as  2a_(P„/m  -  R  )  -  P  or  a.  - 

2o:^.''m,  et.c.  If  we  can  neglect  all  the  two-Gaussian  intogr-?.is 

then  the  Gaussians  become  decoupled.  Consider  now  a  situation  in 

which  <?(t)  varies  in  time  faster  than  ail  other  quantities.  If 

we  analyze  the  behavior  of  the  Eq .  (111.24)  in  the  neighborhood 

of  a  time  t  we  can  expand  <?(t)  =  <>(t  )  -  (3!})(t  ).'3t  )(t-t  ). 

0  0  0  0 

The  exponential  term  e^*^''"'  oscillates  (in  the  neighborhood  of 
t^)  with  the  period  T  =  2  tt/ [  S?  ( t  ^) /3 1  ^)  .  If  this  is  smaller 
than  <.'.e  time  scale  ~  over  which  X.(t),  M(t),  V(t),  a(t)  -and  b(t) 
change  .appreciably  we  can  integrate  the  equation  (III. 24)  from 
t  -  T.'2  to  t  -  T/2  and  obtain  .M(t)  X  (t)  =  V(t).  Thus  the  two 


Gaussian  integrals  which  cause  the  coupling  between  the  Gaussiajis 
disappear  from  the  equation  of  motion  and  the  Gaussians  evolve 
independently.  We  see  that  ps£udo-classical  motion  behaves  just 
like  the  classical  one:  it  tends  to  ignore  forces  that  act  at 
frequencies  vastly  different  than  the  rate  of  change  of  the' 
parameters  being  propagated. 

One  can  derive  an  expression  for  :{>{t)  and  show  that  in  the 
cases  when  the  two  Gaussians  overlap  well 

the  phase  of  the  integrals  is  proportional  to  the  difference 
between  .the  phases  of  the  two  Gaussians.  Within  LHA  these  phases 
are  proportional  to  the  classical  actions  along  the  trajectories 
of  the  centers  of  the  two  Gaussians.  So,  tWo  Gaussians  following 
trajectories  having  classical  actions  that  change  ra'^idly  in 
time,  are  weakly  coupled. 


4G  • 


IV.  The  Choice  of  Initial  Wave  Function 
IV.  1  Introductory.  Remarks 

The  choice  of  the  initial  wave  functlo.n  is  in  principle 
very  simple:  it  must  fit  as  closely  as  possible  the  experimental 
conditions  of  interest.  The  practical  implementation  of  this 
idea  in  the  context  of  QWP  propagation  was  done  in  a  manner  which 
causes  ambiguities  and  (sometimes)  trouble. 

The  first  difficulty  appears  because  of  the  practice  of 
writing  the  initial  wave  function  as  a  sum  of  Gaussians  in  a  way 
that  leaves  us  free  to  choose  certain  parameters  (i.e.  width, 
momentum,  etc.)  almost  at  will.  This  "asymptotic  freedom" 
permits  us  sometimes  to  affect  substantially  and  arbitrarily  the 
final  wave  function;  this  is  not  a  desirable  feature  in  any 
theory . 

The  second  difficulty  is  more  subtle  and  is  common  to  all 
methods  using  a  pre-selected,  basis  set  to  represent  the  wave 
function  throughout  the  collision  process.  A  set  might  be 
flexible  enough  to  represent  the  initial  state  well,  but  be 
incapable  to  describe  the  intermediate  or  the  final  wave  function 
with  the  desired  accuracy.  The  problem  is  particularly 
interesting  in  cases  with  many  channel  final  states  of  the  kind 
that  can  be  intuitively  described  by  multiple  classical 
trajectories  that  cover  different  regions  of  configuration  space. 
Such  situations  cannot  be  represented  by  a  single  Gaussian 
packet.  The  desired  flexibility  can  be  achieved  by  increasing 
the  number  of  Gaussians  used  to  fit  the  wave  function.  There  is 
however  a  limit  t.o  this  and  our  experience,  drawn  from  a  variety 
of  numerical  studies,  .is  that  we  cannot  mindlessly  add  more  and 
more  Gaussians  until  the  results  converge,  since  in  the  course  of 
collision  the  Gaussians  often  overlap  causing  over-completeness; 
when  this  happens  the  differential  equations  propagating  the 
parameters  become  singular  and  intractable. 
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IV. 2  The  Asymptotic  Freedow 

To  understand  how  this  problem  arises  it  is  best  to  exaniine 
several  e.xamples.  The  first  is  the  representation  of  a  planar 
wave  as  a  sum  of  Gaussians,  which  has  been  used  in  all  the  GWP 
diffraction  calculations  published  so  far.^®’^®  We  start  with 
the  identity"*^ 


k  e,xp[i^‘ir]  =  C  j*  dr^  e:<p{(i/Fi)  ^  A  •  (T-'r  ) 

z  ' 

■  •'i*'  \  •1*'  < 

4-  ik-(r-r  )  -  ik-r  ) 


( IV. 1  ) 


where  C  is  a  normalization  constant. 

To  obtain  an  approximate  representation  of  the  planar  wave 
as  a  sum  of  Gaussians  we  discretize  the  integral.  This  gives 


k”^'^  e.xpCi^*^]  =  C  Z  exp{  ( i/h)  [  ('r-'r  )  (?-r  ) 

P  P  P 


k'(r-r  )  -  ik*r  }. 
P  P 


(IV. 2) 


The  number  of  Gaussians  and  their  mean  positions  are  fixed  by  the 
accuracy  we  i.mpose  when  we  represent  the  integral  by  a  sum.  The 
momentum  of  each  packet  is  hk  and  the  phase  is  real  and  given  by 
k'r^.  However  the  method  gives  no  prescription  for  fixing  the 
initial-  values  of  the  width  matrix  A.  It  is  reasonable  to  take 
the  initial  off  diagonal  elements  zero  and  assume  that  the 
diagonal  ones  are  equal,  because  of  the  isotropy  of  space.  These 
decisions  still  leave  th'e  complex  diagonal  element  a  of  the  width 
matrix  unspecified. 


The  e-xisiting  practice  has  been  to  argue  that  since  (IV. 2) 
represents  the  initial  state  well  for  any  reasonable  choice  of  a, 
we  can  use  the  "asymptotic  freedom"  to  select  a  value  of  a  that 
would  make  our  life  simpler.  If  we  plan  to  use  SHM  (which  has 
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been  the  case  so  far)  we  should  select  Imo:  so  that  the  packet 
will  be  narrow  when  it  collides  with  the  hard  wall  of  the 
potential.  This  should  increase  the  accuracy  of  LHA  (which  is 
used  in  SHM).  However  since  Ima(t)  is  controlled  by  the 
equations  of  motion  we^can  not  fix  its  value  at  the  wall  by 
selecting  the  initial  value.  The  practice  has  been  to  use  the 
equation  of  motion  for  Imo£(t)  in  free  space  and  to  select  Ima(O) 
such  that  Iri!a(t)  at  the  wall  location  would  be  large  if  the 
packet  moves  in  free  space.  While  this  gives  some  guidance 
concerning  the  choice  of  Jm£)c(o)  it  leaves  Realo)  unspecified  and 
this  is  taken  to  be  zero. 

1  8 

Unfortunately  detailed  numerical  studies  show  that  this 
choice  of  the  width  does  not  achieve  its  stated  goal:  no  matter 
how  we  choose  Ima(o)  the  potential  broadens  the  packet  beyond  the 
values  for  which  LHA  can  be  safely  applied.  Furthermore,  we  find 
that  the  final  results  depend  sometimes  on  the  choice  of  Ima(o). 
While  in  the  case  of  diffraction  changing  Irta(o)  does  not  lead  to 
large  deviations  from  the  known  quantum  results,  we  feel  rather 
uncomfortable  in  using  such  a  strategy  for  cases  where  the 
"exact"  results  are  unknown. 

Another  example  is  Heller's  integral  representation  of  a 
harmonic  oscillator  wave  function 

T 

J*  dt  exp{-(|^)  (y-y{t))^  +  {  i/h)  p  ( t )  ( y-y  ( t )  ) 

°  (IV. 3) 

+  ( i/2h) (p(t)y(t)  -  p(o)y(o))  +  inut} 

Here  p(t)  and  y(t)  satisfy  the  classical  equations  of  motion  of 
the  oscillator  momentum  and  position,  T  =  2'rT/u,  and  p(o),  y(o) 
are  the  initial  conditions  for  the  momentum  and  position  of  the 
oscillator.  We  can  now  represent  us  a  sum  of  H  Gaussians 

by  discretizing  the  integral.  This  gives 
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«i)n(y)  =  S  exp{-(m<d/2li)[y-y(t  )3^  +  (l/R)  P  ( t^)  [  y-y  ( t^)  ] 
a=  1 

(IV. 4) 

+  (i/2fi)(P(t  )y(t  )  -  p(o)  y(o))  +  Inut  }  , 

C2  C(  OC 

I 

wi  th  t  “  (  2-iTa)  /  {  Ncd )  . 
ct 

The  prescription  tells  us  that  the  points  P ( >  V ( 
classical  trajectory  at  equa.lly  spaced  time  intervals.  The 
initial  p(o),  y(o)  are  hot  however  specified  so  the  phase  of  the 
classical  oscillatory  motion  giving  p{t),  q(t)  is  arbitrary. 

A  similar  situation  occurs  in  the  representation  of  the 
rotational  wave  functions,  where  grotip  theory  tells  us  how  to 
construct  the  wave  function  as  a  sum  of  Gaussians  whose  centers 
are  located  on  the  surface  of  a  sphere.  The  other  parameters  in 
the  Gaussians  remain  at  our ’ disposal .  In  more  general  cases  the 
asymptotic  wave  functions  <j>^(x:0)  are  represented  as  linear 
combinations  of  N  Gaussians 

N 

4)^(x:0)  =  z  c^^.G^(x:{X(o)}^)  (iv.5) 

A=  1 

whose  parameters,  symbolized  in  (IV, 5)  by  {X(o))^,  are  chosen 

before  the  linear  coefficients  C.  are  determined.  The  latter 

An 

are  found  by  minimizing  the  total  energy  of  the  asymptotic 
system;  in  the  course  of  this  minimization  the  parameters 
{X(0)}^  are  frozen.  This  procedure  also  suffers  from  the  fact 
that  it  provides  no  objective  method  for  chosing  {X(0)}^. 

In  what  follows  we  discuss  a  procedure  which  is  more 
efficient  and  more  satisfactory  conceptually  and  practically:  we 
represent  the  initial  wave  function  !j>^(x:0)  by  a  sum  of  Gaussians 
whose  parameters  are  chosen  by  a  non-linear  least  square  fitting 
(LSF)  procedure.  That  is,  we  minimize 
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F{  {X 


A  =  1 


Xd.\  x(x)  i  {({>^(x;  o) 


N 

2  G  (x:{Xl 
A  =  1 


2 


{  IV.6) 


with  respect  to  the  parameters  {X}^.  The  idea  is  so  simple  that 
it  would  not  merit  further  discussion  except  for  the  fact  that  it 
brings  about  a  number  of  dramatic  improvements. 

(a)  The  number  of  Gaussians  required  for  obtaining  a  good 
fit  by  Eq.  (IV.6)  is  much  smaller  than  that  required  by  oth.er 
methods.  Consider  for  example  the  fit  of  a  low  lying  Morse  state 
by  using  Eq.  (IV. 5).  We  can  make  a  reasonable  choice  of 
Gaussians  as  follows.  If  we  assume  that  the  low  lying  Morse 
states  are  nearly  harmonic  we  can  use  Heller's  equation  (IV. 4)  to 
select  the  parameters  (X)^  (i.e.  position,  momentum  and  phase)  in 
G^{x:{X}^).  Taking  linear  combinations  of  these  Gaussians,  like 
in  (IV. C),  we  can  find  the  linear  coefficients  by  minimizing 

the  energy  with  respect  to  them  and  keeping  (X)^  frozen.  We  can 
get  very  good  fits  of  the  low  lying  states  by  usi.ng  eight 
Gaussians.  By  using  Eq.  (IV.6)  we  obtain  an  equally  good  fit 
with  only  two  Gaussians. 

(p)  The  non-linear  least  square  fit  method  has  the 
advantage  that  it  fixes  all  the  parameters  objectively.  The 
number  of  Gaussian  is  predetermined  by  the  choice  of  the  error 
that  we  are  willing  to  tolerate  in  the  initial  wave  function  and 
the  flexibility  required  daring  propagation. 

(y)  It  is  interesting  to  note  that  lowering  the  number  of 
Gaussians  is  not  a  matter  of  efficiency  only.  We  find  that  it  is 
very  difficult  to  propagate  wave  functions  co.m posed  of  a  large 
number  of  coupled  Gaussians  because  In  the  course  of  their 
evolution  they  can  overlap  and  the  set  become. s  overcomplete.  As 
a  result  the  differential  equations  propagating  the  parameters 
become  nearly  singular  and  give  very  large  errors.  As  an  example 
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propagated  a  linear  combination  of  eight  Gaussians  representing 
the  third  Morse  state  and  a  three  Gaussian  non-linear  least 
square  fit  (i.e.  Eq.  IV. 6)  to  the  same  function.  In  the  first 
case  the  computer  fails  to  solve  the  .MEM  differential  equations 
propagating  the  para,meteps,  because  the  matrix  M  coupling  the 
Gaussia.ns  becomes  singular.  The  reason  for  this  seems  to  be  the 
overcompletness  of  the  set,  which  we  detect  by  diagonalizing  the 
matrix  formed  with  the  Gaussian  overlap  integrals:  the  singular 
behavior  of  M  is  always  preceeded  by  the  decrease  of  one  or  more 
of  the  eigenvalues  of  the  overlap  matrix.  A  calculation  using  a 
sum  of  three  Gaussians  to  fit  non-linearly  the  initial  Morse 
state  has  no  difficulty. 

It  is  important  to  note  that  the  non-linear  fitting  is  not 
entirely  free  of  ambiguities,  since  -several  "best  fits"  can  be 
obtained,  depending  on  the  starting  point  and  the  minimization 
strategy  pursued.  Consider  for  e.xample  the  third  Morse  eigen¬ 
function  which  the  non-linear  least  square  fit  program  can 
represent  very  well  by  a  sum  of  three  Gaussians.  Let  us  assume 
now  that  we  decide  to  try  a  four  Gaussian  fit.  We  find  chat  for 
certain  starting  parameters  the  LSF  program  makes  the  amplitude 
of  one  Gaussian  nearly  zero  and  fits  the  wave  function  with  the 
remaining  three  Gaussians.  Even  though  we  get  a  very  good  fit 
this  sum  is  a  very  bad  initial  function  since  the  MEM  program 
cannot  propagate  it;  the  matrix  It  in  Eq  .  (A. 10)  is-  nearly 
singular  because  of  poor  overlap  between  the  nearly  zero 
amplitude  Gaussian  and  the  others.  However  it  is  quite  possible 
to  get  a  satisfactory  four  Gaussian  initial  wave  function  if  we 
constrain  the  width  and  the  normalization  1.11(7'^)  for  each 

Gaussian  to  stay  within  reasonable  limits.  The  MEM  progra.T. 
propagates  this  function  rather  well. 


IV. 3  The  OptlauM  Number  of  Gau^^sians 

In  many  cases  we  would  like  to  use  the  smallest  number  of 
Gaussians,  and  fo;"  this  the  non-llne'ar  LSF  of  the  initial  wave 
function  is  very  helpful.  There  are  however  cases  when  such  a 
choice  would  be  physically  unsound.  Consider  a  Morse  oscillator 
colliding  with  an  atom.  We  can  fit  the  Initial  wave  function 
(the  ground  state  of  the  Morse  oscillator)  well  with  one 
Gaussian.  However,  if  the  kinetic  energy  of  the  incident  atom  is 
comparable,  but  smaller,  than  the  dissociation  energy,  the  final 
state  is  a  linear  combination  of  several  Morse  functions.  One 
Gaussian  cannot  describe  correctly  such  a  wave  fuirction;  it  can 
at  best  give  the  average  energy  transferred  but  we  could  not 
expect  corre'ct  state  occupation  amplitudes.  It  is  therefore  a 
good  idea  to  try  to  fit  the  initial  state  with  several  Gaussians. 
This  is  a  general  situation  in  most  cases  in  which  the  final 
state  is  very  different  from  the  initial  one. 

Another  situation  requiring  a  fit  to  many  Gaussian 
functions  is  that  in  which  the  final  state  has  several  channels 
which  are  qualitatively  different.  One  simple  example  is  atom 
surface  collisions  in  which  surface  trapping  is  of  comparable 
likelihood  with  surface  reflection. 

These  situations  are  too  subtle  and  rich  in  physical 
consequences  to  be  treated  profitably  in  the  general  setting  of 
this  paper.  Several  specialized  studies  of  photodissociation, 
vibrational  excitation  of  diatomics,  curve  crossing  and  surface 
trapping,  which  provide  interesting  and  detailed  illustrations 
for  the  importance  of  chosing  correct  multiple  Gaussians 
representations  of  the  wave  function,  will  be  published  shortly. 
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V.  FINAL  STATE  AXALYSIS 

Generally  scattering  theory  requires  the  koowletige  of  the 

eigenstates  of  the  2ero  order  Hamiltonian  and  the  c,omputation  of 

var-ious  matrix  elements  involving  them..  One  of  the  advantages  of 

the  Gaussian  wave  packet  approach  is  that  we  can  calculate  easily 

the  matrix  elements  needed  for  the  propagation  of  the  wave 
\ 

function:  most  potential  energy  functions  can  be  fitted  to 

polynomials,  exponentials,  Gaussians,  or  to  sums  or  products  of 
such  functions,  so  the  integrals  can  be  done  easily:  the  matrix 
elements  of  the  kinetic  energy  operator  require  the  calculation 
of  integrals  involving  a  .product  of  Gaussians  and  polynomials. 

Part  of  this  advantage  is  however  lost  if  we  must  analyse 
the  scattered  wave  functions  by  calculating  the  matrix  elements 
S  ( X  )'i' ( X  ,  t )  dx  with  the  eigenstates  of  the  final  zero-th  order 
Hamiltonian.  From  a  practical  point  of  view  in  many  situations 
we  don't  have  simple  formulae  for  4>^(x)  and  we  must  generate  them 
numerically,  which  makes  the  calculation  even  more  tedious.  And 
sometimes  are  known  only  very  approximately. 

We  present  below  a  very  simple  and  rather  general  idea  that 
permits  the  analysis  of  the  final  wave  function  by  using  the 
program  that  propagates  the  MEM  equations.  Since  the  overlap  of 
planar  waves  with  Gaussian  functions  is  a  Guassian  in  momentum 
space  we  need  no  special  procedure  for  the  analysis  of  the  final 
translational  state.  We  concentrate  therefore  on  analysing  the 
internal  states  only. 


Let  us  assu.me  that  at  the  time  when  the  projectile 
target  interaction  stopped  we  have  a  scattered  wave  function 
(internal  state)  given  by 


'i'{x:  t^) 


A 


(V.i) 


We  can  use  the  ME.M  equations  to  propagate  this  wave  function  wit.^: 
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the  zero  order  Hamiltonian  of  the  internal  ’states.  If  we  denote 
. . . “  ■■■ 

by  ij){x;t)  the  propagated  wave  function  we  can  easily  show  that 
«  4-i[(i)-iX](t-t'  ) 

C(u)  2  X  dt  e  {x:t)  il)(x:t^)dx  (V.2) 

^0 

satisf ies 


Re  C(u)  =  Z - — 2  ‘  (V..3) 

n  =  0  (w-(j  )  -  X 

n 


and 


(u-u  ) 

Im  C((i))’=  Z  :c  :  - S - ^  .  (V.4) 

n  (u-u  )  e  X*^ 

n 


Here  n  runs  over  the  bound  (internal)  states  of  the  system. 

2 

!C^|  is  the  occupation  of  n-th  state  (i.e.  the  probability  that 

the  scattering  process  takes  the  system  into  its  n-th  state)  and 

is  the  energy  of  that  state.  The  quantity  X  is  at  our 

disposal.  If  we  make  it  much  smaller  than  u  the  peaks  in  ReC(u) 

2 

are  well  separated  and  are  given  by  peak  positions  and 

by  the  peak  height.  The  zeroes  of  Im  C(u)  are  close  to  . 

However,  if  X  is  too  small  then  we  must  propagate  the  wave 

function  il)(x:t^)  for  a  long  time  t  such  that  t  X  >>  1.  A 

compromise  can  be  reached  by  using  an  intermediate  value  for  X 

2 

and  determining  and  by  a  least  square  fit  of  C(cj)  to  the 

forms  (V.3)  and  {V.4). 


From  a  physical  point  of  view  the  quantity  C('J)  is  a 
Green’s  function  for  a  fictitious  absorption  process  (or' 
fluorescence)  which  is  used  to  resolve  the  post-collision  state 
ij/(x;t)  into  spectral  'components. 
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''Sometimes  the  propagation  procedure  is  not 


reliable  e  n  o  u  g  h 


to  be  comfortably  used  for  a  very  long  time.  In  such  occassions. 
we  can  define  C.j,(u)  by  using  Eq.  (V.2)  with  the  upper  integration 
limit  equal  to  T.  It  is  easy  to  show  that 


ReC.j.(c5) 


and 

Im  C.^(u) 


Z  ’C  X[(G-u  ^  X"^e“^^ 

■  n  n 

jy 

[Xcos(u-u  )T  )  sin(u)-(j)  )T]} 

-  n  n  n 


{  V  .  3  ) 


n 

(V.6) 

[Xsin  ) cos  { u-o) } 


We  can  use  these  equations  and  a  least  square  fit  procedure  to 
2 

find  (i)  and  I C  ‘  . 

n  ■  n 

3y  using  a  fast  Fourier  transform  we  find  that  this 
procedure  is  both  efficl'*it  and  reliable.  An  e.tample  was  shown 
in  Fig.  7.  This  was  obtained  as  follows;  we  made  the  linear 
combination 

;;i(.\:t^)  =  a^4)^  -  a^^*^  -  ^  a^^i^ 

where  6  is  the  n-th  .Morse  eigenfunction;  we  then  fitted  tnis 
n 

function  to  four  Gaussians  and  pretended  that  this  is  our  post¬ 
collision  function:  we  propagated  the  Guassians  with  the  Morse 
Hamiltonian  and  the  MEM  equations  (solid  curves).  The  graph 
shows  ReC(u)  calculated  by  using  E.q .  (V.2).  The  least  square  fit 

analysis  of  these  curves  gives  the  eigenvalues  and  spectral 

,2  ,  .2 

composi  tion  ;  a . !  •  • 

0  0  • 


VI  STABILITY  PROBLEMS  IN  THE  PROPAGATION  OF  THE  MEM  EQUATIONS 
V I  •  1  Introductory  re.-narks  ♦ 

In  principle  MEM  would  permit  the  propagation  of  huridreds 
of  Gaussians,  making  many  localized  time  dependent  quantum 
problems  within  the  reach  of  today's  computer  power. 

Unfortunately  our  numerical  experience  has  revealed  some 
limitations  which  are  summarized  in  this  section. 

The  first  limitation  is  a  numerical  instability  in  the 
propagation  of  the  Gaussian’s  widths.  This  was  encountered  by 
Heller  in  his  use  of  SHM  and  he  circumvented  it  by  using  what  we 
call  here  a  P-Z  transformation.  We  show  that,  not  unexpectedly, 
the  same  difficulty  is  present  in  MEM  and  that,  fortunately,  the 
ME.M  equations  can  be  written  in  a  form  which  permits  the 
application  of  the  P-Z  method. 

The  second^  limitation  appears  when  we  attempt  to  use  a 
large  number  of  Gaussians.  We  find  that  in  the  course  of  ti.T.e 
'the  Gaussians  often  evolve  in  a  way  that  makes  one  {or  more)  of 
them  redundant.  When  this  happens  the  .ME.M  equations  become 
singular  and  cannot  be  solved.  Superficially  this  may  seem  a 
pleasant  problem,  to  be  solved  by  reducing  the  initial  number  of 
Guassians.  Unfortunately  the  opti.mu.m  number  of  Gaussians  is  not 
a  uniform  function  of  time:  as  the  collision  proceeds  the  wave 
function  contracts  or  spreads  {in  coordinate  and/or  momentu.-n 
representation)  so  that  the  number  of  necessary  Gaussians  goes  up 
and  down  in  time.  While  at  some  given  time  N  Gaussians  may  be 
too  many  and  cause  trouble,  they  may  be  needed  at  other  times. 

'We  found  no  simple,  general  method  of  dealing  with  this  problem, 
but  we  designed  a  useful  strategy  that  is  present  here. 

Since  in  most  problems  of  interest,  to  us  the  exact  auantum 
solution  is  not  known  we  test  for  errors  in  the  propagat.. 
scheme  by  looking  for  internal  ir. consistencies.  Practically  we 


>=^’n;T,  '• 


‘.-.v-v.. 


sr 


use  three  critJjria:  (1)  we  require  <(|i(t)!ti)(t)>  to  be  time 

indepentleiit  :  (2)  we  reauire  <!i!  (  t )  ;  H  , 'ii  {  t )  >  to  be  time  independent; 

and  (3)  we  require  'S  dt  exp[  +  iut  -  XtJ  <i}i  ( t )  ;  ij)  (  o  )  >  (which 

2 

satisfies  Eqs.  (V.3-4))  to  give  values  of  which  add  up  to 

one  (when  «jj(t)j:]:(t)>  =  1)  and  correct  values  for  the  eigen - 

energies  . 

n 

VI. 2  The  P-Z  Transformation 

The  P-Z  transformation  was  designed  by  Heller  to  solve 
difficulties  connected  with  the  propagation  of  the  width  matri.K 
a.  In  SHM  the  difficulty  appears  in  the  equation  c;  = 

-2oc^/m  -  (1/2)  o^V/SR^  (Eq.  (III. 7))  propa'ga  t  ing  the  width 
parameter  a.  This  has  an  oscillatory  behavior  which  causes  a  lot 
of  trouble  if  we  apply  usual  numerical  methods  (i.e.  Runge-Kutta 
or  pred-ictor-correc  tor )  to  Eq .  (III. 7).  In  the  best  situations 
this  can  be  cured  by  using  an  extremely  small  time  step.  In 
other  cases  erroneous  values  are’ obtained  even  for  tha  smallest 
time  steps.  Our  experience  has  been  that  both  diffraction  and 
curve-crossing  calculations  with  SHM  require  the  use  of  the  P-Z 
method . 

The  ME.M  calculations  carried  out  by  us  so  far  show  that  a 
direct,  numerical  solution  of  the  MEM  equations  lead  to  dramatic 
failures,  'much  more  rapidly  and  frequently  than  in  the  case  of 
SHM ; 


Fortunately  the  P-Z  transform  can  be  applied  directly  to 
the  MEM  equations  if  they  are  written  in  the  proper  form.  Ve 
start  with  the  equations  (A. 10-13)  written  as 

^.1-  — y  1 

.X  =  (  M  )  •  V  H  ?  ( VI  .  1  ) 


The  only  equations  in  the  above  system  that  require  modification 

» 

are  those  containing  the  components  of  X  haviTig  tne  form  s  - 
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2 

2a^/m.  For  example,  in  the  one-aimensional  two-Gaussian  case 
discussed  in  Appendix  A,  there  are  the  components  X_  and 

in  Eq.  (A. 11).  For  three-dimensional  Guassians  these  equations 
have  the  matrix  form 


o: 


A 


(VI  .2) 


where  o:  is  the  three  dimensional  matrix  appearing  in  the  term 

(i/a)  (x  -  R(t))*a  •  (x-R(t))  at  the  exponent  of  each  of  the 

three  dimensional  Gaussian;  F  (t)  is  a  known  function  of  time. 

A 

We  can  remove  the  non-linear  term  2  a^'cs^/m  by  introducing 
two  new  variables  2  and  P  (where  P  is  not  to  be  confused  with  the 
momentum)  through 


The  time  derivative  of  a  is  given  by 

•  •  -1  _ ^-1 

2"  =  "  "  -  "  d(  2  )/dt 

=  ^ 
If  we  now  define 


(VI  .3) 


(VI  .4) 


Z  H  P  /m  .  (VI .5) 

and  use  (VI. 3-5)  in  (VI. 2)  we  obtain 

=  2^  •  (VI  .6) 

The  P-2  methods  solves  (VI. 5)  and  (VI. 6)  and  uses  the 
results  to  compute  a  from  (VI. 3). 


We  find  that  the  use  of  this  procedure  cures  dramatically 
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some  of  the  problems  appearing  when  we  try  to  solve  {VI. 2) 
directly . 

VI . 3  The  Singular  behavior  of  the  MEM  equations 


c 


To  solve  the  MEM  equations  we  must  invert  the  matrix  ^ 
appearing  in  (A. 10).  Its  elements  are  various  .moments  M(3n;Am) 
of  the  Gaussians  used  to  fit  the  wave  function.  For  a  one- 
dimensional  two-Gaussian  wave  function  the  matrix  M  is  given  by 
(A. 12).  The  left-hand  upper  corner  of  that  matrix  is  the  overltip 
matri.x  between  the  Gaussians  used  to  make  the  fit.  This  suggests 
that  if  the  overlap  matrix  becomes  singular  it  may  be  difficult 
to  invert  M.  Empirically  we  find  this  to  be  the  case.  As  we 
solve  the  equations  of  motion  for  the  parameters  we  also  solve 
for  the  eigenvalues  of  the  overlap  matrix.  We  find  that  whenever 
one  eigenvalue  becomes  very  small  the  determinant  of  M  becomes 
small  and  large  propagation  errors  appear.  For  problems  with  a 
small  number  of  parameters  it  is  better  to  diagonal ize 


G.;'  is  minimized.  If  we  do  not  interfere,  t.he 
A  ■  • 


There  are  a  va-riety  of  methods  which  we  use  to  cure  this 

problem.  (a)  In  some  cases  the  problem  is  created  by  the  way  in 

which  the  sum  of  Gaussians  represent  the  initial  state.  As  a 

simple  example  consider  the  case  when  we  want  to  fit  a  very  broad 

Gaussian  G  with  four  narrower  Gaussians  G,  A=l,...4.  We  do  this 

A 

by  varying  the  parameters  in  the  Gaussian  so  that 
E  =  s  ;  G  - 

minimization  program  might  decide  to  vary  the  parameters  in  G, 

and  make  it  identical  to  G,  while-  simultaneously  making  Imy^.... 

Imv,  so  large  that  G.,  G.,  G,  have  very  small  values.  Any 

attempt  to  propagate  the  function  Z  G.  obtained  in  this  way  by 

A  « 

MEM  leads  immediately  to  catastrophic  errors.  One  can  easily 
prevent  the  above  events  by  minimizing  E  and  keeping  Im(y)^  below 
a  preset  valu.e.  The  type  of  behavior  exemplified  above  tends  to 
be  general.  We  find,  for  example,  that  as  we  increase  the  number 
of  Gaussians  used  to  fit  the  second  excited  state  of  a  Morse 


&*-  -  •’ 
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oscillator  an  unconstrained  four  Gaussian  fit  makes  the  amplitude 
of  one  Gaussian  ver'y  small  and  the  propagation  of  the  resulting 
wave  function  fails  very  rapidly.  If  we  constrain  the  amplitudes 
we  get  a  four  Gaussian  fit  which  can  be  propagated  successfully. 

We  note  that  symmetry  can  play  an  equally  important  role. 
For  example,  consider  an  initial  wave  function  which'  is  symmetric 
around  .  A  four  Gaussain  fit  might  use  one  Gaussian  centered 
at  and  two  placed  symmetrically,  and  make  the  amplitude  of  the 
fourth  nearly  equal  to  zero.  If  we  keep  the  centers  of  the 
Gaussians  symmetrical  around  R^  ,  all  four  Gausslans  are  use.d  but 
the  fit  may  be  of  poor  quality  if,  for  example,  the  wave  function 
peaks  at  R^. 

While  such  poor  starts  can  be  easily  identified  and  cured, 
there  are  cases  when  in  the  course  of  its  time  evolution  the 
spatial  extent  of  the  wave  function  shrinks  causing  more  serious 
difficulties.  If  N  Gaussians  are  required  to  fit  the  wave 
function  at  times  when  it  has  a  large  spatial  extent,  they  may  be 
redundant  when  the  function  shrinks.  In  such  a  case  the 
propagation  program  may  either  make  the  Gaussians  linearly 
dependent  or  make  the  amplitude  of  one  of  the  Gaussians  zero.  In 
all  these  cases  we  find  that  the  propagation  gives  large  errors 
or  stops  altogether. 


A  cure  for  this  proble.m  can  be  provided  by  monitoring,  the 
evolution  of  the  overlap  matrix  (or  the  M  matrix)  eigenvalues  in 
time  and  by  removing  one  Gaussian,  when  one  eigenvalue  becomes 
small.  This  can  be  done  by  fitting  the  current  N'  Gaussian  wave 
function  tlj(x;t)  =  ?  6  .  ( x  :  {X(  t ) }  , )  to  N'-l  Gaussians  whose 

<  A  A  A 

parameters  {X  (t)}.  are  fitted  to  minimize 
N  N-1  ' 

6  =  j|  2  G.(x:{X}  )  -  Z  G.(x;{X  ),)i|.  This  removes  the 

^  A  A  n  A  ^ 

problem  but  it  can  create  a  new  one  latter:  if  the  wave  function 
expands  spatially  we  may  find  ourselves  with  insufficient 
Gaussians  to  describe  properly  this  expansion.  We  can  however 
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add'a  Gaussian  by  using  the  reverse  procedure:  fit. the  current 
N’-l  Gaussian  wave  function  having  the  parameters  {X  )  ,  to  a  sum 
of  N  new  Gaussians,  by  adjusting  the  paramters  ^  to  minimize 
S.  We  know  that  addition  of  more  Gaussians  is  needed  when  the 
eigenvalues  of  the  overlap  matrix  are  all  close  to  one. 

This  procedure  is  useful,  but  unfortunately  requires  a 
programmer's  supervision  and  interaction  in  the  course  of 
propagation:  one  cannot  do  research  and  play  tennis 
simultaneously,  and  this  can  only  diminish  the  popularity  of  the 
method. 

VI. 4  The  use  of  frozen  Gaussians 

•  * 

We  should  mention  that  many  of  these  difficulties  are  eased 
by  the  use  of  Gaussian  wave  functions  with  fixed  widths  which 
Heller  calls  frozen  Gaussians.  Since  the  width  is  not  changing, 
the  P-Z  transform  is  not  required.  Empirically  we  also  find  that 
the  numerical  stability  of  the  .ME.M  equation  with  frozen  Gaussians 
is  greater.  At  this  time  our  opinion  is  that  more  complex 
problems  will  be  attacked  more  successfully  by  using  frozen 
Gaussians.  The  lack  of  flexibility  caused  by  the  use  of  a  fixed 
width  can  be  compensated  by  increasing  the  number  of  Gaussians 
(without  necessarily  increasing  the  number  of  equations). 

Acknowledgement 

This  work  is  supported  by  the  National  Science  Foundation 
{ CHE82-06 1 30 )  and  in  part  by  the  Office  of  Naval  Research.  We 
are  grateful  to  Eric  Heller,  Howard  Taylor  and  Walter  Kohn'for 
useful  conversations  regarding  this  work. 


APPENDIX  A 


The  Propagatioa  of  multiple  Gauss  Ian,  wave  functions. 

In  this  Appendix  we  derive  a  set  of  equations  frequently 
used  in  the  text.  They  give  the  MEM  propagation  equations  for 
the  case  when  the  wave  function  is  expressed  as  a  sum  of 
Gauss ians : 


ij)(?{;t)  SZ  G^{x:{a^(t).  y^(t)}.  {R^(t)  .P^(t)})  (A.l) 

A 

with 

G^(x;{a^(t).  y^(t)}  ,{R^{t).  PA(t)})  = 

exp{  (i/h)  [o:^(x-R^(t)  +  P^(x-R^(  t )  )  +  y^]}  •  (A. 2) 


If  we  compare  with  the  general  equation  of  section  II. 3. A  we  have 


•ai 


“a'  ''az 


’'a'  '‘ai 


"a’  ‘•aZ  ■  fA- 


In  order  to  obtain  equations  for  a^,  y^,  R^  and  P^  for  all 
the  Gaussians  we  use  the  equations  (11.25).  The  calculations  are 
lengthy  but  straightforward.  The  results  are  listed  below. 


For  Bp  (in  Eq.  II. 25. a)  corresponding  to  the  width 

parameter  a„  of  the  Gaussian  B  we  have: 
o 


Z  {<(x  -  Rg)^  °bI  2a^/m) 

A 

*<(>‘-«b>^‘’bI°a^<’'a-''a®a  -  ‘’'V*  " 

*«>‘-''b)'“bI<’‘-'*A>'=A>  f2“A''’A^”  -  «A>  * 

+  <(x-Rg)^Ggi  V  G^>}  =  0  .  (A. 3) 

If  Bp  (in  Eq.  II. 25. a)  corresponds  to  the  parameter  y  of 

D 
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trajectories  in  a  different  form.  To  do  this  w,e  use  the  notation 
M{An|Bm)'H  ;dx(x-R^)"  G*  (x-Rg)"'  Gg  (A. 8) 


V(An|Bm)  -  J‘dx(x-R^)"G*(x)  V(x)  (x-Rg)"Gg(x) 
We  can  then  summarize  the  equations  (A3-A7)  as 


(A. 9) 


M 


X  *  V 


where 


X  s 


^A^^"  ■  ^A^A  “  ^  V 

Pg/2m  -  PgRg  -  ihog/m  +  y, 


2«a(Pa/"‘  -  Ha)  -  P, 
2ag(Pg/m  -  Rj,)  +  P 


a.  +  2a*/ra 
A  A' 


“b  ^  2ag/m 


B‘ 


B 


(A. 10) 


(A.n) 


M  5 


M(A0|A0) 

M(A0  1  BO') 

0 

M(A0 1 B1 ) 

M(A0|A2) 

M(  AO 

|B2) 

M(B0| AO) 

M(B0| BO) 

M(BOIAl) 

0 

M(B0! A2) 

M(  BO 

|B2) 

0 

M(A1 |B0) 

M(A1 |A1) 

M(A1 i B1 ) 

0 

M(A1 

|B2) 

M(B1 1  AO) 

0 

M{B1 |A1) 

M{B1 1 B1 ) 

M(B1 i A2) 

0 

M(A2|A0) 

M(A2 1  BO) 

0 

.M(A2|B1) 

M(A2|A2) 

M(A2 

jB2) 

M(B2 1  AO) 

M(B2|B0) 

M(B2| Al) 

0 

.M(B21  A2) 

M(  B2 

|B2) 

V{AO| AO) 

+. 

V{A0 1 Bo) 

V(B0|A0) 

+ 

V(BO| BO) 

V(A1 1  AO) 

+ 

V(A1 1  BO) 

V(B1|A0) 

+ 

V(B1 JBO) 

V(A2 1  AO) 

+ 

V(A2| BO) 

V(B2 |A0) 

+ 

V(B2|  BO)/ 

Ai\ 


Another  useful  form  Is 


M  •  X  =  V 

•S#  Am 


with 


^1  ^  ^(«a) 


Xg  ^  V(Rg) 


X3  -  3V(R^)/3R^ 
X4  +  3V(Rgy3Rg 


Xg  +  (1/2)3^V/3R^ 
Xg  t  <l/2)3^V/3Rg 


and  ^  is  obtained  from  ^  by  replacing  V(An|Bm)  with 

V(An(Bm)  =  J-dx{(,x-R^)"  G*{x)[V(x)  -  V(Rg)  - 
-(3V(Rg)/3Rg)(x-Rg)  >  ( 1 / 2 ) ( 3 ^ V ( R ^ ) / 3 R ^ ^ " R ^ ^ ^ 
•(x-Rb)”  Ggix)} 

Thus  SHM  is  obtained  by  taking  "v  =  o 
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Figure  1  . 


Figure  2. 


Figure  3  , 


Figure  4  . 


Figure  Captions 

The  "length",  defined  as  l(t)  =•  (n/ 2 1  mo:  ( t )  )  ^ ^  ,  of  a 
normali-zed,  initially  narrow,  low  energy  Gaussian 
wave  packet  propagated  in  a  Morse  potential  {Eq. 
(.III. 12a))  using  the  MEM  equations  (Eq.  (III. 5)) 

( - ).  l(t)  is  in  units  of  X~"  and  time  is  in 

1  /  9  _  1  _  1 

units  of  T  =  (2,11/0)  '  TT  X  ,  where  X  =  1.846  A 

2  —28  2 

and  D  =  4.334  a.ra.u.  A  / 10  sec  are  the  .Morse 

potential  range  parameter  and  well  depth, 

respectively,  and  m  =  0.5  a .  m .  u .  is  the  reduced  .mass 

of  the  oscillator.  For  values  of  1  higher  than  .22. 

the  LHA  equation  for  the  force  is  in  error  by  more 

than  lOSj.  Initial  values  of  the  wave  packet 

- 1 4 

parameters  are:  Re(c«)  =  0.0  a.m.u./lO  sec:  Im(a) 

- 1 4  - 1  4 

=  8.0  a.m.u./lO  sec;  P  =•  0.0  a.m.u.  A/10  sec: 

R  =  0.20  A:  Re(y)  =  0.1878  a.m.u.  A^/IO"^'^  sec:  and 

Im(r)  =  -0.0696  a.m.u.  A^/IO"^"^  sec. 

The  force  (in  units  of  DX)  e.'terted  on  the  center  of 
the  wave  packet  whose  parameters  are  defined  In  Fig. 
1,  propagated  using  the  MEM  equations  (Ill.Sb) 

( - ‘)  and  the  LHA  equations  (111.7b)  ■ 

( - )  . 

The  Morse  potential  averaged  over  Gaussians  of 
different  width,  <G  v;g>.'<G  G>  (in  units  of  D),  as 
a’function  of  R  (in  units  of  X  M  using  i  =  0.60 

( - --),  1  =  0.329  (  -^ - ),  and  1  =  0.147 

( - )  :  we  also  plot  the  .Morse  potential  V(R)  { - ). 

The  potential  energy  V(R(t))  (in  units  of  D)  as  a 
function  of  time  for  the  wave  packet  of  Fig.  1 

propagated  using  the  .MEM  equations  ( - )  and  the 

LHA  equations  (--  -  "“)• 


Figure  5.  The  quantum  energy,  <G|H|G>  (in  units  of  D),  as  a 
function  of  time  for  the  wave  packet  of  Fig.  1 
propagated  using  the  MEM  equations  ( — and  the  LHA 
equations  (--  — 

Figure  6.  The  "classical  energy,"  <G1V|.G>  +  P^/2m  (for  MEM)  or 

V(R)  +  P^/2m  (for  LHA),  as  a  function  of  time,  for  the 
*  wave  packet  of  Fig.  1.  propagated  using  the  MEM 

equations  ( - )  and  the  LHA  equations  (-- - t  --). 

Figure  7.  The  real  part  of  the  Fourier  transform  of 

<i(i  (  X  ;  t )  1 1(1  ( X  :  t^ )  >  ,  Eq  .  (V.2),  where  the  initial  wave 
function  ili(x:t^)  =  0.5  (<j>^{x)  +  4>j(x)  +  <1>2(^)  ^  4*3(5')) 
is  represented  by  four  wave  packets,  (})^(x)  are  the 
Morse  eigenstates,  and  ij)(x;t)  is  propagated  using  the 

MEM  equations  ( - ),  and  the  uncoupled  IGA  equations 

{--  -  --).  The  peaks  of  ReC((i))  are  related  to  the 

probability  of  being  in  eigenstate  4>j^{5')  (V.3). 

In  this  plot  T  »  5.0,  X  =  1.06. 

Figure  8.  The  square  of  the  projection  of  the  wave  function  (Eq. 

III. 25)  propagated  by  IGA  onto  the  Morse  eigenstates, 

2 

i-e-i  I  1  4*  ( t )  >  I  •  versus  time  for’  n  =  0  ( - ), 

n=l  (--),  n=2  (-  •  -)  and  n=3  (--  -  --). 
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